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Preface 


After  being  selected  for  the  doctoral  program  at  AFIT,  my  first  task  was  to  choose  a 
dissertation  topic  and  a  research  advisor.  Before  I  even  arrived  in  Ohio  I  had  already  asked 
Dr.  Bridgman  if  he  would  take  me  on  as  his  student,  in  addition  to  his  many  duties  as 
Associate  Dean  for  Research.  He  agreed  to  lead  me  in  a  couple  of  special  studies  in  his 
research  specialty:  nuclear  fallout  modeling.  After  my  coursework  and  exams  were  com¬ 
pleted,  we  agreed  that  nuclear  cloud  rise  modeling  was  an  area  we  both  wanted  to  learn 
more  about. 

I  would  like  to  thank  Dr.  Bridgman  for  taking  the  time  to  guide  me  in  the  area  of 
nuclear  fallout  modeling.  His  skill  as  a  teacher  was  known  to  me  from  my  time  as  a  mas¬ 
ter's  student  a  few  years  earlier.  His  reputation  as  an  excellent  researcher  became  evident 
as  I  began  speaking  with  others  in  the  nuclear  weapons  effects  community.  His  search  for 
unbiased  truth  was  driven  home  in  our  weekly  meetings.  His  lifetime  commitment  to  edu¬ 
cation  was  seen  in  his  commitment  to  my  success  and  those  around  me.  Finally,  his  ability 
to  joke  with  me  while  keeping  me  motivated  toward  the  goal  of  graduating  was  summed 
up  by  his  favorite  line  to  me:  “I  got  my  degree!’’ 

I  would  also  like  to  thank  Dr.  Oxley  and  Lt.  Col.  Hartley  whom  I  enjoyed  as  teachers 
as  well  as  committee  members.  I  am  also  honored  for  having  Dr.  John  serve  as  my  Dean’s 
Representative  for  my  dissertation.  I  am  also  thankful  for  the  other  “nukes”  in  the  depart¬ 
ment,  both  students  and  teachers,  which  helped  me  out  during  my  years  at  AFIT. 

The  research  using  DELFIC  and  especially  the  particle  rise  modeling  would  not  have 
been  accomplished  without  the  advice  and  feedback  of  Dr.  McGahan  at  SAIC.  I  greatly 
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appreciate  the  many  hours  that  he  and  Dr.  Cockayne  spent  with  me,  both  listening  and  teaching.  I 
was  very  fortunate  to  have  access  to  Dr.  Cockayne  who  has,  in  my  experience,  the  best  command 
of  both  theory  and  experiments  related  to  nuclear  weapons  effects.  I  would  like  to  thank  LTC 
Byers  of  DNA  who  arranged  access  to  the  TASS  results  and  provided  computer  time  for  plotting 
the  results.  I  am  grateful  to  Tom  Dunn  and  Dr.  Bacon  for  helping  me  understand  the  megabytes  of 
data  involved  and  the  workings  of  the  TASS  code. 

My  wife  Carol  and  my  three  beautiful  daughters,  Jenna,  Danielle,  and  Christina,  deserve  a  lot 
of  the  credit  for  my  success.  Their  encouraging  words  and  ability  to  listen  to  my  troubles  helped 
me  get  through  some  rough  periods  without  losing  hope.  I  thank  them  for  their  sacrifice  of  family 
time  and  a  “normal”  job  schedule,  especially  near  the  end  of  my  research.  Most  of  all  I  thank  my 
Creator  God  who  gave  me  the  desire  and  ability  to  do  the  work.  Out  of  all  my  studies,  nothing 
compares  to  the  fulfillment  and  worth  of  what  I’ve  learned  over  the  past  few  years  about  Him  and 
his  love. 

Vincent  J.  Jodoin 
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Abstract 

Nuclear  cloud  rise  and  growth  to  stabilization  is  investigated.  The  Cloud  Rise  Module 
of  the  Defense  Land  Fallout  Interpretive  Code  (DELFIC)  is  reviewed  and  improved.  Both 
the  modeling  of  bubble  rise  dynamics  and  particle  rise  dynamics  are  addressed.  Improve¬ 
ments  are  made  in  predicting  the  stabilized  particulate  cloud  for  nuclear  fallout  purposes. 
The  results  are  compared  to  US  atmospheric  nuclear  weapons  tests  and  to  the  output  of 
MAZ/TASS  hydrocode  simulations. 

The  improvements  to  the  bubble  rise  dynamics  include  using  a  single  term  entrain¬ 
ment  equation  and  constants  for  the  entrainment  and  eddy  viscous  drag  parameters.  The 
validation  showed  a  slight  improvement  was  obtained  with  the  improved  model  in  predict¬ 
ing  the  cloud  top  values  for  54  US  near-surface  bursts. 

The  improvement  to  the  particle  rise  model  consisted  of  replacing  the  1-D  particle  rise 
with  a  2-D  analytical  flow  field  model.  The  model  was  based  on  Hill’s  spherical  vortex. 
This  improvement  provided  a  radial  distribution  to  the  radioactive  dust  cloud.  The  larger 
fallout  particles  (>200  microns)  achieved  a  higher  altitude  due  to  simulating  the  circula¬ 
tion  within  the  cloud.  Dose  rate  contours  for  some  cases  will  be  altered  due  to  this  change 
in  vertical  distribution  from  the  1-D  model. 
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NUCLEAR  CLOUD  RISE  AND  GROWTH 


I.  Introduction 

Background 

Nuclear  cloud  rise  and  growth  is  typically  modeled  in  two  parts:  1)  nuclear  cloud  rise 
and  growth  to  “stabilization”  and  2)  diffusive  growth  and  transport  in  the  atmosphere  after 
stabilization.  The  second  part  requires  a  definition  of  the  stabilized  cloud  from  the  first 
part. 

This  report  describes  efforts  to  improve  the  cloud  rise  and  growth  to  stabilization  pre¬ 
diction  capabilities  of  bubble,  or  single-cell,  cloud  rise  models.  A  single-cell  model 
assumes  constant  properties  throughout  the  cloud  (e.g.  temperature)  as  opposed  to  hydro¬ 
codes  which  use  a  multicell  approach  that  produces  spatially  varying  properties.  An  exam¬ 
ple  of  bubble  or  cloud  rise  models  is  the  Cloud  Rise  Module  (CRM)  of  the  Defense  Land 
Fallout  Interpretive  Code  (DELFIC).  The  Defense  Nuclear  Agency  (DNA)  is  currently 
updating  DELFIC,  the  Department  of  Defense’s  (DoD’s)  primary  research  reference  fall¬ 
out  code,  which  was  last  published  in  1979.  Being  a  research  code,  DELFIC  contains  as 
much  physics  as  practical  while  reducing  the  number  of  empiricisms.  As  a  reference  code, 
DELFIC  is  used  by  many  other  operational  codes  wliich  either  mimic  or  directly  use 
DELFIC  results  for  systems  studies.  Science  Applications  International  Corporation 
(SAIC),  in  particular  Dr.  Joseph  McGahan,  is  directing  the  update  of  DELFIC.  The  focus 
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of  this  study  is  only  on  improving  the  CRM  of  DELFIC,  and  it  does  not  address  improve¬ 
ments  in  the  cloud  transport  model. 

Brief  Physical  Description.  When  a  nuclear  weapon  detonates  near  the  earth’s  sur¬ 
face,  the  fireball  contacts  the  soil  on  the  ground  and  generates  a  large  cloud  of  dust  which 
rises  and  grows  in  the  atmosphere.  After  an  initial  ballistic  phase,  the  fireball  begins  to  rise 
buoyantly  like  a  bubble  in  water.  The  incompressible  phase  starts  once  the  pressure  inside 
the  cloud  comes  close  to  equilibrium  with  ambient  atmospheric  pressure.  The  upward 
motion  of  the  cloud  distorts  the  initial  spherical  shape  and  sets  up  a  toroidal  vortex  motion 
within  the  cloud  (See  Figure  1  on  page  3).  The  cloud  entrains  ambient  air  and  water  vapor 
as  it  rises.  The  transfer  of  potential  energy  within  the  cloud  to  kinetic  energy  of  rise,  fol¬ 
lowed  by  transfer  to  turbulence,  influences  the  rise  and  growth  of  the  cloud,  and  ultimately 
the  location  of  the  dust  within  the  cloud.  Modeling  the  distribution  of  the  dust  within  the 
cloud  is  important  because  it  contains  the  residual  radioactivity  of  the  weapon  “fallout”. 
The  position  of  the  dust  is  a  function  of  time  during  both  cloud  rise  and  transport  to  the 
ground  and  it  determines  the  radiation  hazard  on  the  ground. 

Motivation  for  the  Research.  The  CRM’s  definition  of  the  stabilized  radioactive  dust 
cloud  in  DELFIC  is  problematic  for  two  reasons.  The  first  reason  is  the  errors  that  exist  in 
the  CRM’s  bubble  rise  equations.  These  errors  include  both  theoretical  and  programming 
errors.  The  second  reason  is  the  simplistic  way  in  which  the  CRM  models  the  particulate 
cloud  in  one  dimension.  This  modeling  may  be  the  reason  for  the  discrepancy  between 
DELFIC  results  and  the  data  from  hydrocodes. 

Huebsch’s  Recommendations  on  Bubble  Rise  Modeling.  Huebsch  was  the  origi¬ 
nal  author  of  the  CRM  (22).  Huebsch’s  1975  critique  of  Normen  ’s  1970  version  of  the 
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Figure  1:  Toroidal  circulation  within  the  radioactive  cloud  from  a  nuclear  explosion 

(13:29) 


CRM  stated  that  the  bubble  rise  model  produced  physically  unrealistic  results  or  results 
that  disagreed  with  the  experimental  data  (24:iii).  Upon  closer  examination  of  the  model, 
Huebsch  proposed  the  errors  were  due  to  the  formulation  of  the  equations  for  the  cloud’s 
mass,  velocity,  temperature,  and  dimensions.  He  found  that  some  of  the  terms  in  the  equa¬ 
tions  contained  theoretical  errors. 
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Huebsch  presented  an  amended  set  of  bubble  rise  equations  which  removed  both  the 
theoretical  errors  and  output  discrepancies.  He  made  further  recommendations  on  how  to 
improve  the  CRM  of  DELFIC.  He  also  presented  one  way  in  which  stabilized  cloud  mod¬ 
els  could  be  categorized,  described  later  in  this  report  (See  Table  1  on  page  11). 

In  response  to  Huebsch ’s  critique,  Norment  published  a  validation  study  of  a  modified 
CRM  (46).  Norment  did  change  the  momentum  (or  velocity)  equation,  as  Huebsch  had 
recommended,  but  not  the  entrainment  equation.  He  adjusted  some  of  his  empirical 
parameters  as  another  method  of  reaching  better  agreement  with  experimental  data.  He 
pointed  out  which  parameters,  using  his  set  of  equations,  most  influenced  changes  in  the 
results.  Huebsch  recommended  keeping  only  a  single  term  of  the  entrainment  equation 
based  on  current  theory.  Noiment  argued  against  such  a  move,  desiring  rather  to  have 
results  matching  test  data  than  equations  matching  current  theory.  A  model  which  does  not 
violate  current  theory  and  matches  the  test  data  is  needed  in  research  level  codes. 

McGahan’s  Recommendations  on  Particle  Rise  Modeling.  McGahan  suggested 
using  vortex  theory  and  a  two-dimensional  representation  of  the  particles’  spatial  distribu¬ 
tion  to  replace  the  limited  1-D  modeling  in  the  CRM  (34).  Vortex  theory  could  be  used  to 
provide  a  two-dimensional  flow-field  in  the  rising  cloud.  By  tracing  the  flow  of  particles 
during  the  cloud’s  rise,  the  position,  velocity,  and  acceleration  of  particles  could  be 
tracked  up  to  stabilization.  While  experimental  results  to  validate  particle  rise  history  are 
not  available,  results  of  a  two-dimensional  (many  cell)  hydrodynamics  code  are  available 
for  comparison. 

A  goal  of  this  research  is  to  refine  the  vertical  distribution  of  all  particles,  especially 
the  larger  particles  modeled  in  DELFIC.  Recent  studies  have  shown  that  DELFIC  cur- 
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rently  stabilizes  tbe  larger  particles  at  much  lower  altitudes  than  the  hydrocodes  predict 
(33).  Vortex  theory  combined  with  bubble  rise  theory  would  be  a  better  method  of  calcu¬ 
lating  the  stabilized  radioactive  dust  cloud  in  the  research  level  code  used  by  the  nuclear 
effects  community.  Accounting  for  the  spatially  and  temporally  varying  flow  field  (vortex 
theory)  will  provide  a  better  way  of  tracking  the  particles.  This  approach  will  predict  a 
vertical  variation  which  should  be  different  than  the  current  prediction  of  DELFIC’s 
CRM.  In  addition,  it  will  definitely  provide  a  radial  distribution  of  particles  that  is  not  cur¬ 
rently  available  in  DELFIC’s  CRM.  Not  only  will  this  tracking  give  a  better  definition  of 
the  particle  distribution  with  time,  but  also  give  a  dynamic  representation  of  the  particles. 

ftObkE 

The  problem  statement  is  simply:  To  determine  the  best  available  theory  for  cloud  rise 
and  growth  using  a  single  bubble  model,  and  to  analyze  the  effect  of  using  vortex  theory 
for  particle  distribution  within  the  rising  cloud. 

Scope 

This  research  includes  two  main  sections:  bubble  rise  modeling  and  particle  rise  mod¬ 
eling.  The  bubble  rise  modeling  focuses  on  correcting  the  current  versions  of  DELFIC’s 
CRM,  improving  the  model  based  on  Huebsch’s  recommendations,  and  validating  the  cor¬ 
rected  and  improved  versions’  results  to  atmospheric  nuclear  test  observations.  The  sec¬ 
ond  section,  particle-rise  modeling,  requires  a  more  substantial  change  to  the  procedure  in 
which  the  positions  of  nuclear  cloud  particulates  are  computed.  A  comparison  can  then  be 
made  to  current  DELFIC  procedures  for  computing  particle  positions  and  also  to  the  out¬ 
put  of  a  hydrodynamics  code. 
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The  CRM  corrections  were  straight  forward  and  consisted  of  checking  the  subroutines 
line-by-line  for  consistency  with  the  code’s  documentation.  Each  line  of  the  CRM  was 
reviewed  and  the  logic  flow  checked.  The  improvements  consisted  of  reverting  to  a  single¬ 
term  entrainment  equation  and  changing  two  parameters  to  constants  (the  entrainment  and 
eddy  viscous  drag  coefficients).  The  best  values  for  these  constants  were  then  determined. 

The  improvements  and  corrections  were  validated  against  experimental  measurements 
using  the  observed  cloud  tops  of  US  atmospheric  nuclear  weapons  tests.  The  observed 
cloud  bases  were  also  examined.  The  comparison  presented  is  a  relative  (or  fractional) 
error  of  calculation  to  observation.  A  set  of  54  US  atmospheric  nuclear  tests  was  used  for 
the  comparison.  The  yields  (in  equivalent  tons  of  TNT)  ranged  from  7.8  tons  to  15  mega¬ 
tons  (MT).  The  shots  in  the  validation  set  included  a  wide  range  for  the  height  of  burst 
(HOB).  They  included  a  slightly  buried  burst,  a  surface  burst,  some  near  surface  bursts, 
and  some  pure  air  bursts.  For  all  cases,  actual  measured  values  of  atmospheric  properties 
(temperature,  pressure,  humidity,  and  winds)  at  the  time  of  the  shot  were  inputted  to  the 
codes. 

The  improvements  in  the  particle  rise  modeling  of  the  radioactive  dust  cloud  required 
entirely  new  methods  for  tracking  trace  particles  in  two  dimensions.  An  initial  2-D  grid  of 
trace  particles  were  used  to  describe  the  initially  uniform  distribution  of  mass  in  the  cloud. 
This  replaced  a  1-D  stack  of  disks  which  were  loaded  uniformly  in  the  vertical  dimension. 
Each  trace  particle  was  tracked  up  through  the  rise  history  of  the  cloud  with  both  vertical 
and  radial  velocity  components  to  determine  a  2-D  distribution  of  mass  in  the  stabilized 
cloud.  This  replaced  the  rise  of  disks  with  only  a  vertical  component  of  velocity.  A  study 
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was  made  of  how  the  initial  soil  loading  is  distributed  within  the  cloud  before  it  starts  its 
rise.  Also  noted  was  the  effect  of  the  vortex  flow  with  particle  size. 

Particle  positions  within  the  cloud  have  not  been  experimentally  measured  but  have 
been  calculated  by  hydrocodes.  Contour  plots  of  the  cloud’s  density  were  chosen  as  the 
best  method  for  comparing  two-dimensional  descriptions  of  the  radioactive  dust  cloud 
with  two-dimensional  hydrocode  output.  When  comparing  the  stabilized  particulate  cloud 
to  current  one-dimensional  CRM  output,  the  vertical  distribution  of  the  mass  was  chosen 
for  comparison.  A  study  of  the  radial  distribution  of  particles  in  the  cloud  using  two- 
dimensional  flow  in  DELFIC  was  performed  and  implications  of  the  distribution  on  fall¬ 
out  contours  is  discussed. 

The  test  cases  chosen  for  particle  rise  modeling  included  the  DELFIC  1979  documen¬ 
tation  test  case  (48:34-50)  and  a  range  of  cases  from  the  Single  Burst  Nuclear  Cloud  Data¬ 
base  generated  by  the  hydrocodes  (9: 15).  These  cases  ranged  in  scaled  heights  of  bursts 
(SHOB  =  HOB/yield1^)  from  0  to  120  fr/kt1^  and  yields  from  20  to  800  kilotons. 

Assumptions 

As  with  any  design  study,  this  work  required  certain  assumptions  and  limitations.  It 
was  accepted  that  certain  studies  previously  performed  in  determining  some  initial  condi¬ 
tions  and  model  parameters  of  the  CRM  need  not  be  repeated  in  this  study.  Since  the 
changes  made  to  the  entrainment  equation  were  substantial,  the  model  parameters  directly 
associated  with  it  were  varied.  Although  some  initial  conditions  and  other  model  parame¬ 
ters  were  investigated,  none  were  altered. 


7 


The  test  data  used  for  comparison  (cloud  tops)  were  known  more  precisely  than  other 
values  for  US  atmospheric  tests,  such  as  cloud  bases  and  radii.  No  assumption  is  made  as 
to  the  error  in  the  observed  cloud  top  values.  Likewise,  the  actual  atmospheric  soundings 
(temperature,  pressure,  humidity,  and  wind;  obtained  were  assumed  to  be  the  best  input 
for  the  test  shots  even  though  some  error  exists  in  their  values  because  of  instrumentation 
(1;10;39),  timing,  and  location  of  the  measurements. 

One  may  not  assume  that  the  hydrocode  data  represents  the  radioactive  cloud.  The 
hydrocodes  do  not  distinguish  between  radioactive  and  nonradioactive  dust.  DELFIC,  on 
the  other  hand,  is  solely  concerned  with  radioactive  dust.  Care  must  then  be  taken  when 
comparing  DELFIC’s  particulate  cloud  to  the  hydrocode  results  because  of  this  distinc¬ 
tion.  The  hydrocode  data  used  for  comparison  was  not  produced  as  pan  of  this  study,  but 
rather  is  pan  of  a  database  maintained  by  DNA  for  those  in  the  nuclear  effects  community. 

General  Approach 

The  sequence  of  major  steps  accomplished  as  part  of  this  research  will  now  be 
reviewed  in  chronological  order.  To  begin,  the  1979  version  of  DELFIC  was  modeled  in 
Mathematics  (63),  a  higher-level  programming  language.  The  theoretical/empirical  equa¬ 
tions  of  the  CRM  were  programmed  in  Mathematics  and  solved  numerically.  This  gave  a 
baseline  from  which  to  depart.  Once  the  equations  were  programmed  identically  to  the 
stated  theory  of  DELFIC,  the  results  were  compared  to  DELFIC’s  FORTRAN  output. 
This  identified  FORTRAN  programming  errors  that  were  known  to  exist  (38:53-57)  and 
uncovered  others  that  had  yet  to  be  found.  It  also  showed,  qualitatively,  how  those  pro¬ 
gramming  errors  affected  the  output  of  the  CRM. 
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The  next  major  step  was  to  analyze  the  recommendations  of  Huebsch  that  were  dis¬ 
carded  by  Norment.  The  purpose  of  this  step  was  to  correct  the  equations  being  used  so 
that  they  no  longer  violated  current  theory.  The  theory  that  best  reduced  the  number  of 
empiricisms,  yet  still  reproduced  the  test  data,  was  determined.  A  check  was  then  made 
against  the  original  test  data  previously  used  to  validate  the  CRM.  Once  this  phase  was 
accomplished  the  first  haif  of  the  project  was  complete. 

The  work  on  the  two-dimensional  particle  rise  modeling  started  with  a  review  of  vor¬ 
tex  theory.  In  particular,  it  included  theory  about  the  flow  in  and  around  a  rising  bubble. 
The  theoretical  equations  were  coded  in  Mathematics  for  symbolic,  numeric,  and  visual 
manipulation  and  for  a  better  physical  understanding.  Once  the  theory  was  understood,  a 
review  of  its  application  to  dust  cloud  modeling  was  undertaken.  Next,  additional  FOR¬ 
TRAN  subroutines  were  written  to  implement  the  theory  into  DELFIC’s  CRM.  While 
Mathematica  is  a  very  useful  tool  to  benchmark  a  program  written  in  FORTRAN  or  simi¬ 
lar  languages,  the  cost  for  these  advantages  is  longer  execution  time.  It  is  currently  too 
slow  for  use  in  an  operational  code  on  today’s  computers. 

The  new  method  of  particle  rise  was  then  used  to  produce  a  database  similar  to  that  of 
the  hydrocode  output.  The  two  sets  of  results  were  then  compared  using  contour  plots  of 
the  dust  densities  at  various  times.  The  vortex  modified  code  was  then  compared  to  the 
one-dimensional  particle  rise  of  the  current  CRM.  The  comparison  was  made  of  the  verti¬ 
cal  distribution  of  mass  in  the  stabilized  cloud  since  the  current  CRM  does  not  define  a 
radial  distribution.  Finally,  the  radial  distribution  of  the  two-dimensional  model  was 
reviewed  and  implications  of  the  vertical  distribution  on  fallout  contours  were  found. 
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Sequence  of  Presentation 


The  research  is  now  presented  in  more  detail.  Chapter  II  begins  with  an  overview  of 
the  bubble  rise  theory  in  DELFIC.  It  details  the  equations  which  were  corrected, 
improved,  and  validated  to  give  the  rise  characteristics  of  the  nuclear  cloud  using  an 
entraining  bubble  assumption.  Chapter  HI  then  reviews  how  the  rise  characteristics  can  be 
used  with  vortex  theory  to  solve  the  flow  field  in  and  about  the  cloud.  This  flow  field  is 
used  to  trace  the  history  of  particles  which  end  up  defining  the  stabilized  nuclear  particu¬ 
late  cloud.  Chapter  IV  examines  the  results  of  both  the  bubble-rise  and  the  particle-rise 
test  runs.  The  implications  of  these  results  are  discussed  in  Chapter  V. 
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H  Bubble  Rise  Dynamics 


The  purpose  of  this  chapter  is  to  explain  the  theory  of  bubble  rise  dynamics  and  state 
the  work  done  by  the  author  in  this  general  area.  The  chapter  is  divided  into  three  main 
sections.  The  first  section  is  a  short  overview  of  one  way  in  which  cloud  rise  models  can 
be  classified.  The  next  section  is  much  larger  and  explains  the  level  of  effort  used  in  this 
work.  The  final  section  is  a  breakdown  of  the  steps  taken  to  correct,  improve,  and  then  val¬ 
idate  this  part  of  the  cloud  rise  model. 

Levels  of  Effort 

There  are  three  main  levels  of  effort  for  a  cloud  rise  model,  as  detailed  by  Huebsch 
(24: 18-19)  (See  Table  1  on  page  1 1).  A  first  level  effort  is  one  that  obtains  the  stabilized 
cloud  results  through  empirical  fits.  It  develops  time  histories  of  cloud  height,  dimensions, 
and  temperature  as  functions  of  yield  and  height  of  burst.  The  effort  is  not  responsive  to 
atmospheric  inputs.  The  empirical  fits  could  be  derived  from  atmospheric  test  observa¬ 
tions  or  from  the  output  of  a  second  or  third  level  effort.  The  K-division  DNA  Fallout 
Code  (KDFOC)  is  an  example  of  fallout  codes  which  use  this  approach  to  model  the  stabi¬ 
lized  cloud  (15). 


Table  1 :  Levels  of  Effort 

Level  Variables 

1  Height,  dimensions,  and  temperature 

2  Above  plus  turbulence,  mass,  and  partitions  of  mass 

3  Above  plus  pressure  and  density 


Method 
Empirical  fits 
One  zone  cloud 
Multiple  cell  cloud 
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A  second  level  effort  is  one  that  uses  a  bubble,  or  single  zone,  cloud  model  assumption 
to  solve  conservation  equations  for  the  stabilized  cloud  properties.  The  bubble  contains 
global  properties  for  the  entire  cloud  (e.g.  a  single  temperature  everywhere  in  the  cloud). 
This  effort  calculates  the  time  histories  of  the  cloud’s  height,  dimensions,  temperature,  tur¬ 
bulent  kinetic  energy,  and  fractions  of  mass  in  the  form  of  air,  water  vapor,  condensed 
water,  and  soil.  The  model  has  yield  and  height  of  burst  (HOB)  as  inputs.  It  is  also  respon¬ 
sive  to  atmospheric  conditions.  That  is,  the  cloud  rise  depends  on  the  pressure,  tempera¬ 
ture,  water  content,  and  other  variables  derived  from  them,  of  the  atmosphere  through 
which  the  cloud  is  rising.  This  level  of  effort  uses  the  continuity  equations  to  obtain  a  cou¬ 
pled  set  of  nonlinear  ordinary  differential  equations,  with  time  as  the  only  independent 
variable.  This  set  is  then  solved  to  give  the  time  histories  mentioned. 

A  third  level  effort  is  a  full  hydrodynamics  code  solution  to  the  cloud  variables.  Local 
horizontal  and  vertical  variations  within  the  cloud  are  modeled  in  a  multicell  representa¬ 
tion  of  the  cloud.  This  could  be  a  2-D  or  3-D  level  hydrocode  which  also  provides  infor¬ 
mation  such  as  pressure  and  density  of  the  individual  cells.  The  independent  variables 
include  both  time  and  position  coordinates.  Ambient  atmospheric  properties  are,  of 
course,  included  in  the  calculation. 

Existing  hydrocodes  include  MAZ  and  TASS.  These  are  typically  run  on  supercom¬ 
puters.  The  MAZ/TASS  tandem  solves  the  cloud  dynamic  equations  for  the  compressible 
(<1  min.)  and  incompressible  (1  min.  to  1  hr.)  phases  of  the  cloud’s  history,  respectively. 
These  codes  have  additional  parameters  which  must  be  determined/chosen,  which  are 
unique  to  this  higher  level  of  modeling. 
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The  DoD  research  reference  fallout  code,  DELFIC,  contains  three  main  modules:  the 
Cloud  Rise  Module  (CRM),  the  Diffusive  Transport  Module  (DTM),  and  the  Output  Pro¬ 
cessor  Module  (OPM)  (See  Figure  2  on  page  14).  The  CRM  has  two  main  functions.  The 
first  is  to  calculate  the  solution  of  the  bubble  rise  history  of  a  single  zone  cloud,  as 
described  previously  for  a  second  level  effort.  The  second  function  is  to  simulate  the  parti¬ 
cle  rise  to  describe  the  stabilized  particle  cloud,  which  is  needed  as  input  for  the  DTM. 
This  second  function  will  be  described  in  more  detail  in  “Particle  Rise  Dynamics"  on 
page  40. 

History  of  DF.T .FTC’s  Cloud  Rise  Module.  The  original  transport  models  of  the  1950s 
used  an  empirically  computed  stabilized  cloud  as  input  for  their  fallout  calculations 
(61: 11;  18).  These  empirical  fits  provided  the  dimensions  and  height  of  the  cloud  once  it 
stopped  rising.  An  assumed  distribution  of  radioactive  mass  within  the  cloud  was  then 
made  that  produced  adequate  agreement  with  measured  fallout  on  the  ground. 

In  the  late  1950s,  it  was  agreed  by  the  nuclear  effects  community  that  using  physics  in 
modeling  the  cloud  rise  and  growth  (instead  of  empiricism)  was  the  major  area  missing  in 
the  current  models.  In  the  early  1960s  Huebsch  developed  a  model  for  surface  bursts  over 
water  for  the  Naval  Radiological  Defense  Laboratory  (NRDL)  that  included  a  physical 
modeling  of  the  cloud  rise  and  particle  g  owth,  instead  of  using  an  empirical  stabilized 
cloud  (20).  He  then  modified  his  model  to  include  a  soil  burden  for  land  surface  bursts  and 
implemented  it  in  the  original  version  of  the  DoD  research  reference  fallout  code  DELFIC 
(22).  Huebsch  ’$  model  was  the  original  CRM  of  DELFIC  and  was  published  as  Volume  HI 
of  seven  volumes  in  1967  (22;41;42;43;55;59;62). 
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Figure  2:  Top  Level  Flow  Chart  of  DELF1C 


The  CRM  has  seen  many  revisions  since  1967  (See  Table  2  on  page  15).  In  1970,  Nor- 
ment  (45)  revised  the  CRM  and  among  other  changes,  completely  removed  the  particle 
growth  capability.  In  1975  Huebsch  (24)  critiqued  the  1970  CRM-Revised  of  Norment. 
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Huebsch  showed  that  the  1970  equations  violated  the  conservation  of  energy.  Huebsch’s 
recommendations  were  addressed  in  Norment’s  1977  validation  of  CRM-Revised  (46). 


Table  2:  History  of  DELFIC’s  Cloud  Rise  Module 


Author 

Date 

Title 

Huebsch 

1964 

Development  of  a  Water-Surface-Burst  Fallout  Model:  The  Rise  and 
Expansion  of  the  Atomic  Cloud. 

Huebscn 

1967 

The  Department  of  Defense  Land  Fallout  Prediction  System.  Vol¬ 
ume  m.  Cloud  Rise. 

Norment 

1970 

The  Department  of  Defense  Land  Fallout  Prediction  System.  Vol¬ 
ume  ID.  Cloud  Rise.  Revised. 

Huebsch 

1975 

Analysis  and  Revision  of  the  Cloud  Rise  Module  of  the  Department 
of  Defense  Land  Fallout  Prediction  System  (DELFIC). 

Norment 

1977 

Validation  and  Refinement  of  the  DELFIC  Cloud  Rise  Module. 

Norment 

1979 

DELFIC:  Department  of  Defense  Fallout  Prediction  System.  Volume 

I  -  Fundamentals. 

In  the  1977  version  of  the  CRM,  Norment  revised  the  equations  and  model  assump¬ 
tions  so  they  no  longer  violated  the  conservation  of  energy.  One  change  occurring  in  this 
revision  was  the  removal  of  the  virtual  mass  factor  on  the  momentum  (velocity)  equation. 
The  effect  of  removing  this  factor  is  negligible  for  the  times  of  interest  in  DELFIC.  How¬ 
ever,  the  virtual  mass  factor  is  needed  when  considering  bubbles  starting  from  rest  to  limit 
their  acceleration  to  twice  that  of  gravity  (32:4). 

Finally,  in  1979,  Norment  published  the  last  well-documented  revision  of  DELFIC 
(47;48).  Other  modifications  to  DELFIC,  including  conversion  from  FORTRAN  IV  to 
FORTRAN  77,  have  been  accomplished  since  then  (31).  These  other  modifications,  how- 
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ever,  have  not  been  well  documented  or  produce  questionable  results,  and  therefore  the 
1979  version  was  selected  as  the  baseline  for  this  effort  (34). 

Initial  Conditions.  The  CRM’s  computations  are  started  at  some  time  other  than  the 
time  of  detonation.  In  the  CRM,  initial  conditions  are  calculated  as  a  set  of  empirical  fits  to 
nuclear  test  data  (functions  of  yield  and  HOB).  This  provides  the  initial  size,  mass,  height, 
temperature,  and  velocity  at  the  initial  time  for  the  cloud  rise  calculations.  The  selection  of 
the  initial  time,  approximately  1  to  20  seconds  after  detonation,  occurs  before  the  majority 
of  the  cloud’s  rise,  but  after  pressure  equilibrium  with  the  atmosphere  has  been  achieved. 
This  allows  an  incompressible  model  assumption.  By  the  initialization  time,  the  cloud  has 
cooled  sufficiently  to  assume  radiative  losses  of  energy  are  negligible.  The  initial  time  also 
is  needed  for  the  entraining  bubble  assumption,  since  the  very  early-time  rise  of  the  fire¬ 
ball  lacks  entrainment  of  ambient  air  (20;  14). 

A  thorough  review  of  the  initial  conditions  was  completed  to  determine  if  any  should 
be  altered  in  this  study.  Most  were  found  to  be  based  on  observation,  with  only  a  few 
based  on  conjecture.  The  two  initial  conditions  that  contained  the  most  uncertainty  con¬ 
cerning  their  origin  were  the  initial  mass  and  temperature  of  the  soil  in  the  cloud.  Over  the 
many  revisions  of  DELFIC,  however,  the  value  for  the  soil’s  temperature  has  been 
adjusted  to  ensure  reasonableness  in  its  value  for  all  yields  (28).  The  initial  soil  mass  has 
not  been  revised  since  the  original  DELFIC,  but  was  found  to  differ  only  slightly  from  the 
more  modem  fit  by  Carpenter  (5).  The  difference  was  minor  when  compared  to  the  CRM 
assumption  that  the  cloud  obta;ns  immediate  loading  of  the  entire  soil  mass  at  initial  time. 
Therefore,  none  of  the  initial  conditions  were  changed  from  the  way  they  were  calculated 
in  the  1979  version  of  DELFIC. 
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Ordinary  Differential  Equations.  Starting  with  the  initial  conditions,  the  CRM  uses 
the  conservation  equations  of  mass,  momentum,  turbulent  kinetic  energy,  and  heat  to  solve 
for  eight  basic  variables  (see  Table  3  on  page  17).  The  derivations  of  these  equations  are 
contained  in  Huebsch’s  1964  document,  and  revised  by  Norment  in  his  later  publications 
(45;46).  The  equations  are  transformed  into  eight  first-order  nonlinear  ordinary  differen¬ 
tial  equations  (ODEs)  in  the  basic  variables.  These  variables  are  vertical  velocity,  height, 
temperature,  turbulent  kinetic  energy,  mass,  and  fraction  of  mass  in  the  soil,  vapor,  and 
condensed  water  phases.  Given  these  variables,  all  other  variables  of  interest  can  be  deter¬ 
mined.  All  variables  solved  for  are  global,  in  that  the  cloud  model  is  a  single  bubble  with 
uniform  properties  throughout. 


Table  3:  Basic  Variables  Solved  for  in  the  Cloud  Rise  Module  ODEs 
Symbol  Variable  Balance  Equation 


u  vertical  velocity  of  rise  momentum 

z  vertical  height  of  the  cloud  center 

T  temperature  of  the  cloud  heat 

E  turbulent  kinetic  energy  density  turbulent  kinetic  energy 

m  total  mass  of  the  cloud  (air,  soil,  water,  &  vapor)  mass 

s  ratio  of  condensed  soil  in  cloud  mass  to  that  of  air 

x  ratio  of  water  vapor  in  cloud  mass  to  that  of  air 

w  ratio  of  condensed  water  in  cloud  mass  to  that  of  air 


A  short  derivation  for  some  of  the  ODEs  follows.  This  derivation  assumes  that  no  con¬ 
densed  water  is  present  in  the  cloud.  A  more  complete  derivation  may  be  found  in  the 
works  of  Huebsch  (20;21;22). 
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Velocity  Equation.  The  bubble’s  velocity  equation  is  derived  by  equating  its  rate 
of  change  of  momentum  to  the  sum  of  its  buoyancy  force  and  its  eddy  viscous  drag  force 


,  rate  of  change ,  r  buoyancy ,  ,  eddy  viscous , 

{  }={}+{  } 

of  momentum  force  drag  force 


(1) 


w  2/t,vu 

-(mu)  =  V{pe-p)g - m 


CPe 


(2) 


where 


v  =  Max  (u,  J2E)  is  the  characteristic  velocity  (m/s) 
m  =  total  cloud  mass  (air,  vapor,  water,  and  soil)  (kg) 
u  =  rate  of  cloud  rise  (m/s) 

V  =  cloud  volume  =  m/p  (m3) 
pc  =  ambient  air  density  (kg/m3) 
p  =  cloud  density  (kg/m3) 

g  =  acceleration  due  to  gravity  (9.8  m/s2) 

Jl'2  =  eddy  viscous  drag  parameter 
E  =  turbulent  kinetic  energy  density  (J/kg) 

Hc  =  cloud  vertical  thickness  (m) 
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Then  using  the  gas  law  for  air  and  water  vapor  one  can  rewrite  the  buoyancy  term  using 
the  following  relation 


(3) 


where 

T*  =  r1  is  cloud  virtual  temperature  (K) 
l+*/£ 

Te  =  Te—~  is  ambient  virtual  temperature  (K) 

1  ^ 

|3'  =  ^  +  w  is  ratio  of  cloud  gas  density  to  total  cloud  density 

T  =  cloud  temperature  (K) 

Tc  =  ambient  temperature  (K) 

x  =  cloud  water  vapor  mixing  ratio  (ratio  of  water  vapor  mass  to  dry  air  mass) 
xt  =  ambient  water  vapor  mixing  ratio 

E,  =  18/29  is  ratio  of  the  molecular  weight  of  water  vapor  to  that  of  dry  air 
w  =  cloud  condensed  water  mixing  ratio 

Then  one  can  rearrange  terms  to  give  the  differential  equation  for  the  time  rate  of  change 
of  cloud  velocity.  The  first  term  on  the  right  in  Eqn  (4)  is  the  buoyancy  term  while  the  sec¬ 
ond  term  includes  the  eddy  viscous  and  entrainment  drag  effects 
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(4) 


du 

dt 


(  7* 

-^P’  -  1 

T  * 

t 


\ 


8- 


(  2  1  dm 

HrT*P  mdt 
c  1  e 


Energy  Equation.  The  time  rate  of  change  of  turbulent  kinetic  energy  per  unit 
mass,  E,  is  obtained  by  balancing  the  sources  and  sinks  of  kinetic  energy  in  the  cloud  dur¬ 
ing  a  small  time  interval  dt 


^  source  by  eddy  + 
viscous  drag 


source  by  inelastic-^ 
collision 
entrainment 


diluted  by  ^  ^  dissipated  ^  ^ 

entrainment  to  heat 


—  -  'tir  7*  u2ldm 

dt  =  2T*P  h.  2  m  dt 


ent 


c  1  dm 
- Em T, 


-£ 


(6) 


ent 


where 


—  =  rate  of  change  of  mass  due  to  entrainment  (kg/s) 

1  tent 

k*  (2£) 3/2 

e  s  - — - is  turbulent  kinetic  energy  dissipation  rate  per  unit  mass  (J/kg  s) 

Hc 

k$  =  a  constant 

Temperature  Equation.  The  temperature  equation  can  be  obtained  by  the  conser¬ 
vation  of  energy(heat)  during  a  time  interval  dt. 
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beat  absorbed  in 


heat  gained  from 


heat  used  to  warm^  +  ^  - - - }  -  {  - * -  }  (7) 

entrained  air  dm  adiabatic  expansion  dp  turbulent  energy  dissipation 


'  T  ' 

l^jhdf 

KT.  J 


dm  +  m 


[ ?p(T)dT - 


RaT*p4EJ  s  £mpdt 


(8) 


where 

Cpa(7)=  specific  heat  of  entrained  air  (dry  air)  (J/kg  K) 

cp(7>=  mass  weighted  mean  of  the  specific  heats  of  air  and  water  vapor  (J/kg  K) 

Ra  =  gas  constant  of  air  (287  J/kg  K) 

p  =  pressure  at  cloud  center  altitude  (Pa) 

Then  using  the  hydrostatic  equation 


(9) 


one  can  transform  the  expansion  term,  divide  through  by  dt,  and  rearrange  terms  to  get  the 
time  rate  or  change  of  cloud  temperature.  The  first  term  inside  the  brackets  on  the  right 
below  accounts  for  adiabatic  expansion,  the  second  accounts  for  entrained  air,  while  the 
third  is  for  turbulent  kinetic  energy  dissipation  to  heat. 
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(10) 


dT 

dt 


r 

(  T  \ 

p' 

lTe  j 

1  dm 

?p(n 

T  +&u  + 

1 1 

(3 'm  dt 
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Mass  Equation.  The  equation  for  the  rate  of  change  of  mass  due  to  entrainment  as 
derived  by  Norment  is  shown  below. 


where 

S  =  surface  area  of  the  cloud  (m2) 

H  =  entrainment  parameter 

This  version  of  the  entrainment  equation  is  much  different  than  the  original  version  and 
that  suggested  by  Huebsch  in  Eqn  (12) 


dm 

dt 


ent 


S„ 

=  m-nv 


(12) 


Norment ’s  Eqn  (11)  has  two  extra  terms  inside  the  parentheses  on  the  right  representing 
temperature  effects  and  hydrostatic  expansion.  These  terms  and  the  leading  multiplier  on 
the  right  arise  due  to  Norment’s  assumptions  in  his  derivation  (45).  Some  of  Norment’s 
assumptions  are  no  longer  consistent  with  the  current  CRM,  and  Eqn  (12)  is  therefore  the 
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preferred  form  of  the  equation.  Huebsch  has  stated  that  either  form  of  the  entrainment 
equation  can  be  used  without  violating  the  other  conservation  equations. 


Drv  Equations.  Hie  above  four  equations  are  then  combined  with  four  equations 
that  define  the  change  in  mixing  ratios  of  water  vapor,  condensed  water,  and  soil,  and  the 
definition  of  velocity.  The  entire  set  of  eight  ODEs  used  in  the  current  CRM  are  listed  as 
Eqns  (4),  (6),  (10),  (11),  (13),  (14),  (15),  and  (16). 


dx  1  +  x  +  s ,  .  I  dm 


(13) 


where  5  denotes  the  cloud  soil  mixing  ratio  (ratio  of  soil  mass  to  dry  air  mass). 


dw 

Hi 


1  dm 
dt 


ent 


dx 

dt 


(14) 


ds  _ 

11+jt  1  dm 

1  +X  +  S+  w( 

) 

dt  = 

P'  1  +  xe  m  dt 

ent 

m  l 

.  S  +  HV 

(15) 


where  F  denotes  the  rate  of  soil  mass  leaving  the  cloud  (kg/s). 


dz 

dt 


=  u 


(16) 


where  z  denotes  the  cloud  center  height  (m). 
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Wet  Equations.  Some  of  the  equations  listed  above  have  to  be  altered  to  include 
the  effects  of  latent  heat  Since  the  equations  were  derived  from  conservation  relation¬ 
ships,  extra  terms  are  needed  for  the  source  of  energy  released  when  the  water  vapor  in  the 
cloud  condenses.  This  occurs  when  the  saturation  temperature  for  the  pressure  of  the  cloud 
has  been  reached.  The  release  of  latent  heat  can  be  significant  in  some  cases  (20: 10) 
(“Appendix  B:  Latent  Heat"  on  page  129).  The  derivation  of  this  set  of  equations  can  be 
found  elsewhere,  and  it  suffices  to  list  just  those  equations  that  need  altering  (20;21;22). 
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cpRaT2j 


' 
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~guT*( ,  ,  Lx_\  e.1  g“ 
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(17) 


where  L  denotes  the  latent  heat  of  vaporization  of  water  or  ice  (J/kg). 
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(19) 


™=  (l+x/^)-^^  f  S- 

xdt  P  T^dt  v  ^  R  T 


Solution  of  the  QDE's.  The  model  uses  the  set  of  eight  ODEs  and  solves  them 
numerically  using  a  fourth-order  Runge-Kutta-Gili  method  (4:363).  By  using  3  time 
domains  from  initial  time  to  stabilization  time  (f,  to  t,  +  1  sec,  r,+  1  sec  to  r,  +  100  sec,  and 

tj+ 100  sec  to  tst)  with  3  different  time  steps  (1/32, 1/4,  and  2  1/2  sec),  the  CRM  solves  the 
set  of  ODEs  for  each  phase  of  the  solution.  A  somewhat  different  set  of  variables  (base, 
top,  radius,  base-rate,  top-rate,  and  radiai-rate,  T,  and  p)  is  calculated,  printed,  and  stored 
in  an  array  at  various  times  during  the  cloud  rise  to  fully  define  the  history  of  rise  (See 
Table  4  on  page  26).  The  set  of  variables  stored  is  used  to  run  the  later  subroutines  that 
distribute  the  particles  in  the  vertical  dimension.  Plots  of  the  basic  variables,  as  generated 
from  a  typical  run  are  shown  (See  Figure  3  on  page  27). 

Model  Parameters.  In  the  above  equations  there  exist  two  parameters  that  were 
altered  in  this  work  to  be  consistent  with  the  literature  and  to  best  match  observed  cloud 
tops:  parameters  for  entrainment  (|i)  and  eddy  viscous  drag  (ki).  These  parameters  are 

both  currently  functions  of  yield,  W,  (See  Figure  4  on  page  28).  The  original  version  of  the 
CRM  had  these  two  parameters  as  constants:  p=0.2  and  ^=0. 1 .  These  were  both  changed 
by  Norment  to  become  functions  of  yield  in  his  1970  revision.  He  then  changed  them 
again  in  his  1977  validation  to  the  functions  plotted  in  Figure  4  on  page  28.  The  increase 
in  the  slope  of  |i  past  15  MT  is  questionable  since  there  were  no  atmospheric  tests  above 
this  yield  to  validate  against. 
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Figure  4:  CRM  Parameters  for  Entrainment  fyi)  and  Eddy  Viscous  Drag  (k2)  as  a 
Function  of  the  Log  of  the  Yield  (W  in  kt) 


Modeling  Fallout  During  Cloud  Rise.  The  change  in  cloud  mass  during  a  time  step  Jl 
must  also  account  for  the  loss  of  soil  mass,  F,  that  may  settle  out  during  the  rise. 
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The  settling  rates  of  each  particle  size  group  must  be  calculated  at  each  time  step  to  com-  j 

I 

pute  their  contribution  to  the  value  of  F.  A  significant  fraction  of  the  bubble  rise  computa¬ 
tion  is  taken  up  in  computing  this  value.  The  contribution  to  the  change  of  mass  due  to  soil 
loss,  when  compared  to  the  change  due  to  entrained  air,  was  found  to  be  negligible  by  the 
author.  All  non-airbursts  in  the  validation  set  were  used  to  compare  the  effect  of  soil  fall¬ 
out  on  the  basic  variables.  The  effect  of  ignoring  the  soil  loss  was  hardly  noticeable  on  any 
burst  and  was  well  below  a  1%  change  for  all  variables.  The  loss  of  condensed  water  was 
also  modeled  in  the  original  CRM,  but  was  considered  negligible,  and  was  no  longer  cal¬ 
culated  in  Norment’s  CRM-Revised. 

Effect  of  Wind  Shear.  Huebsch  has  shown  that  the  effect  of  wind  shear  on  cloud  rise 
is  important.  The  CRM  models  wind  shear  by  a  replacement  of  variables  in  the  entrain-  j 

ment  equation.  A  physical  explanation  of  the  second  term  on  the  right  below  may  be 
found  in  Huebsch ’s  work  (23:2-5).  Basically  the  effect  of  shear  on  the  cloud  is  to  increase 
the  surface  to  volume  ratio,  thereby  increasing  entrainment. 

(2I) 

where 

k6  =  a  constant 

Rc  =  horizontal  radius  of  the  cloud  (m) 

vs  =  wind  velocity  difference  between  the  top  and  base  of  the  cloud  (m/s) 
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Cloud  Shape  and  Volume.  The  volume  of  the  cloud  can  be  determined  at  any  time  by 
the  ideal  gas  law  for  air,  given  the  temperature,  pressure,  and  mass  of  the  cloud.  Its  dimen¬ 
sions  are  then  determined  by  an  assumption  of  its  shape. 

The  original  CRM  assumed  a  sphere  until  the  cloud  top  reached  a  user-defined  tropo- 
pause  altitude.  The  cloud’s  vertical  radius  was  then  frozen  and  any  further  growth  was 
accommodated  by  horizontal  expansion.  Norment’s  CRM-Revised  took  a  different 
approach.  The  vertical  thickness  of  the  cloud  was  assumed  to  be  linearly  dependent  on  the 
cloud’s  rise  distance.  The  radius  was  thi  n  determined  assuming  a  spheroidal  shape. 
Though  Norment  claimed  this  would  produce  an  oblate  cloud,  Huebsch’s  1975  critique 
showed  that  the  cloud  could  become  prolate  at  late  times. 

Norment’s  response  to  Huebsch’s  critique  was  to  fix  the  cloud  shape  as  an  oblate 
spheroid  of  eccentricity  0.75  based  solely  on  observation  and  calculated  from  the  cloud 
volume.  Once  the  cloud  stopped  rising,  however,  the  cloud  would  stop  growing  vertically 
and  expand  only  radially,  becoming  more  oblate.  This  is  the  current  way  the  CRM  com¬ 
putes  the  cloud’s  dimensions.  Although  this  method  seems  to  perform  better  than  the  pre¬ 
vious  methods,  it  still  does  not  allow  for  the  “collapse”  of  some  clouds  where  the  base 
continues  to  rise  after  the  top  has  stopped  rising.  It  also  does  not  allow  the  top  to  overshoot 
and  then  oscillate  about  a  final  stabilized  altitude.  Both  of  these  phenomena  have  been 
observed  in  US  nuclear  test  events. 

Termination  of  Cloud  Rise  and  Expansion.  Once  the  cloud  has  stopped  rising,  the 
CRM  checks  two  conditions  to  determine  horizontal  stabilization  (47:26-27).  When  either 
of  the  two  conditions  are  met,  the  cloud  expansion  is  terminated  and  the  definition  of  the 
particle  cloud  is  calculated  for  input  to  the  DTM.  The  first  condition  is  that  the  radial 
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growth  rate  has  slowed  down  below  a  predetermined  value  which  is  a  function  of  yield. 
The  other  is  that  the  turbulent  kinetic  energy  density  in  the  cloud  is  reduced  to  ambient 
levels  (also  a  function  of  yield)  for  the  height  of  the  cloud. 

Work  Done  on  DBLFTC’s  Bubble  Rise  Subroutines 

With  the  above  explanation  of  the  current  version  of  the  CRM,  the  work  done  by  the 
author  on  the  bubble  rise  modeling  will  now  be  presented.  This  section  of  the  chapter  is 
broken  down  into  three  subsections.  The  first  subsection  describes  the  corrections  needed 
in  the  CRM  subroutines  to  match  the  stated  theory.  The  second  subsection  details  the 
improvements  implemented  by  the  author  in  the  model.  The  final  subsection  reviews  the 
validation  effort  used  to  check  the  author’s  corrections  and  improvements. 

A  complete  listing  of  the  DELFIC  code  used  in  this  research  is  on  file  in  the  Depart¬ 
ment  of  Engineering  Physics,  Air  Force  Institute  of  Technology,  Wright-Patterson  AFB, 
OH,  45433.  The  code  includes  the  original  1979  version  obtained  by  the  author  along  with 
the  corrections  and  improvements  described  in  this  report.  All  changes  to  the  1979  version 
are  distinguished  by  comments  within  the  code  containing  the  author’s  initials  ‘vjj’. 

Corrections  Needed  in  the  1979  DELFIC  CRM.  In  recovering  the  1979  version  of 
DELFIC  from  the  documentation,  McGaban  had  to  add  a  MAIN  program,  since  the  docu¬ 
mentation  was  lacking  one.  The  CRM  then  compiled  with  no  changes  on  a  386-SX-16 
personal  computer  using  Microsoft  FORTRAN  5.1.  To  compile  on  a  SPARC  station  2 
using  SUN  FORTRAN  1.4,  two  separate  FORMAT  statements  had  to  be  edited.  These 
minor  changes  allowed  the  1979  version  of  the  CRM  to  be  compiled  without  any  warnings 
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or  errors.  Similar  changes  were  needed  to  compile  the  entire  DELFIC  program  including 
the  DTM  and  OPM. 

The  1979  documentation  contains  a  test  case’s  input  and  output  which  were  used  to 
verify  that  the  code  was  running  correctly  (48:34-50).  After  the  CRM  was  compiled  suc¬ 
cessfully,  die  program ’s  output  was  verified  with  that  listed  in  Volume  II  of  the  1979  docu¬ 
mentation.  Additional  post-processors  were  built  by  the  author  with  Mathematics  to 
display  the  CRM  variables  as  a  function  of  time.  DELFIC  is  fast  on  today’s  computers 
with  the  entire  test  case  taking  only  20  seconds  on  a  486-33  personal  computer! 

In  studying  the  1979  version,  a  few  errors  were  found  both  in  the  document  pertaining 
to  the  theory  and  in  the  source  code.  A  discussion  of  these  discrepancies  and  the  correc¬ 
tions  needed  to  remedy  them  are  documented  in  a  previous  report  by  this  author  (25).  One 
correction  needed  in  the  source  code  so  it  would  match  the  1979  documentation  concerned 
the  entrainment  equation.  This  change  made  a  noticeable  difference  in  the  output  of  the 
test  case.  While  the  remainder  of  the  corrections  didn’t  change  the  test  case  output  signifi¬ 
cantly,  they  may  affect  other  scenarios  to  a  greater  extent. 

The  discrepancy  mentioned  in  section  2.3.5  (concerning  a  switch  on  the  specific  heat 
calculation  of  the  cloud),  on  page  11  of  this  author’s  previous  report  (25),  has  now  been 
resolved.  The  switch  is  needed  in  the  calculation  due  to  the  derivation  of  the  temperature 
equation.  The  soil  temperature  is  at  a  different  temperature  than  the  gaseous  cloud.  It  is 
assumed  to  stay  at  this  lower  temperature  until  the  gaseous  cloud  cools  to  this  tempera¬ 
ture,  at  which  time  both  soil  and  gases  are  cooled  at  the  same  rate.  The  interested  reader 
should  refer  to  the  original  CRM  document  (22). 
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Improvements.  The  single  term  form  of  the  entrainment  equation  (See  Eqn  (12))  has 
been  regarded  as  the  correct  form  in  the  literature  since  Taylor  (58),  and  is  still  accepted 
by  the  community  today  (32;54).  Norment  in  his  revision  to  the  original  CRM,  attempted 
to  derive  the  entrainment  equation  from  the  ideal  gas  law  aud  empirical  fits  to  observed 
cloud  behavior.  Huebsch  critiqued  Nonnent’s  version  of  the  entrainment  equation,  consid¬ 
ering  it  invalid  and  the  cause  of  output  discrepancies  (24: 12).  His  recommendation  was  to 
return  to  the  single  term  version  of  the  equation.  Huebsch  also  asserted  that  the  entrain¬ 
ment  (|i)  and  eddy  viscous  drag  (k2)  parameters  should  be  constants,  not  parameters  of 
yield  as  Norment  had  redefined  them. 

Upon  investigating  Nonnent’s  derivation,  Huebsch ’s  recommendations,  and  Nor- 
ment’s  subsequent  changes,  it  was  decided  that  Huebsch ’s  form  of  the  entrainment  equa¬ 
tion  was  best.  When  Norment  derived  his  version  of  the  equation,  he  used  some  rigid 
empirical  fits  to  the  cloud’s  vertical  thickness  as  a  function  of  rise  distance.  He  basically 
abandoned  these  empiricisms  in  his  changes  to  the  CRM’s  approach  to  cloud  shape  in 
1977,  but  chose  not  to  alter  his  entrainment  equation.  This  makes  his  derivation  question¬ 
able  at  best,  and  the  single  term  form  (as  it  still  appears  in  the  literature)  should  be  used. 

The  literature  also  supports  Huebsch  in  stating  that  the  entrainment  parameter,  by  def¬ 
inition,  is  a  constant  (12:7).  There  are  similar  physical  arguments  that  support  Huebsch’s 
recommendation  that  the  eddy  viscous  drag  parameter  too  should  be  a  constant  (20: 18- 
19).  In  his  critique,  Huebsch  mentioned  that  Norment’s  combination  of  increasing  jJ.  and 
lowering  k2  with  yield  compensate  somewhat  for  each  other.  This  explains  why  agreement 
with  observation  is  still  achieved.  It  will  be  shown  later  that  better  agreement  to  atmo- 


33 


spheric  nuclear  test  observations  can  be  achieved  using  the  single  term  entrainment  equa¬ 
tion  and  setting  |i  and  k2  to  constants. 

The  reason  Norment  may  have  achieved  acceptable  results  with  his  yield  dependent 
fits  for  the  parameters  was  due  to  their  compensating  effects  on  cloud  rise.  As  the  reader 
can  see  from  Figure  4  on  page  28,  )i  is  raised  while  k2  is  lowered  as  yield  increases  past 
one  kiloton.  The  result  is  to  increase  the  effect  of  entrainment  while  lowering  the  effect  of 
eddy  viscous  drag  on  the  equations.  When  Norment  reviewed  Huebsch’s  recommenda¬ 
tions,  he  tried  a  single  term  equation  by  setting  |i  and  k2  to  0.25  and  0.1  respectively. 
Although  Huebsch  stated  these  values  as  his  best  guess,  he  did  say  some  testing  was  prob¬ 
ably  required  to  determine  their  best  values.  Norment  did  not  report  varying  these  values 
when  testing  the  single  term  entrainmer*  equation.  The  result  of  using  the  above  values 
was  an  over  prediction  of  entrainment  and  a  poor  match  to  observed  cloud  tops. 

Upon  further  investigation,  the  discrepancy  may  be  found  in  the  cloud  shape.  Hueb¬ 
sch’s  original  model  used  a  spherical  shape.  Norment’s  1977  revision  was  using  a  spheroi¬ 
dal  shape.  This  implies  a  larger  surface  to  volume  ratio  for  Norment.  Therefore,  to 
maintain  the  same  proportionality  factor  m  front  of  the  velocity  term  in  the  entrainment 
equation,  one  must  reduce  (I  for  a  spheroid.  In  fact,  the  best  value  found  here  for  |i  was 
0.12,  which  will  be  discussed  in  “Results  and  Comparisons”  on  page  72. 

Validation.  The  effect  of  correcting  and  improving  the  cloud  rise  equations  are  deter¬ 
mined  by  running  actual  atmospheric  nuclear  test  cases.  The  results  are  then  validated 
against  recorded  test  data  from  these  events.  Input  data  files  were  created  using  the  best 
available  unclassified  sources  which  describe  the  test  events.  Similarly,  observation  data 
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was  gathered  from  published  unclassified  reports.  Finally,  the  figure  of  merit  used  by  Nor- 
ment  was  used  in  this  research  to  judge  the  performance  of  the  changes. 

Norment’s  Validation  Set.  In  his  1970  revision  of  the  CRM,  Norment  included  an 
appendix  which  contained  tabulations  and  plots  from  simulations  of  US  atmospheric 
events  (45: 155-166).  Although  he  used  56  shots,  he  only  identifies  the  54  unclassified 
ones.  These  54  shots  are  the  set  used  in  this  research  for  validation  without  excluding  the 
lowest  yield  shot  as  Norment  did  in  his  1977  validation  of  the  CRM  (See  Table  5  on 
page  36  and  Figure  5  on  page  38)  (46: 18).  Since  Norment  states  that  DELFIC  is  valid  for 
shots  from  1  ton  to  100  MT  (47;48),  this  7.8  ton  event  is  included  in  this  author’s  valida¬ 
tion  studies.  These  54  events  include  39  airbursts,  13  near  surface  bursts,  1  surface  burst, 
and  1  slightly  buried  burst.  The  yields  included  14  sub-kiloton  shots  and  3  megaton  range 
shots.  Only  5  of  the  shots  were  from  the  Pacific  Proving  Grounds,  including  the  four  high¬ 
est  yield  shots.  The  highest  yield  of  the  Nevada  Test  Site  shots  was  74  kilotons. 

Cloud  Top  Observation  Data.  Since  Norment  did  not  cite  any  written  reports  as 
the  source  of  inputs  or  observation  values  for  his  validation,  it  is  impossible  to  know  the 
reliability  of  his  data.  DAS  A  1251-EX-I  and  -II  were  used  in  this  research  for  both  the 
event  data  (yield,  surface  height,  and  height  of  burst)  and  the  cloud-top  observation  data 
(16;  17).  This  choice  provides  consistency  between  the  different  series  (e.g.  Castle  and 
Plumbbob).  It  also  lists  the  modem  best  estimates  of  the  values  and  provides  a  basis  for 
repeatability  for  future  research  done  in  this  area.  Most  event  data  values  agreed  with 
those  used  by  Norment  in  1970. 

Weather  Incut  Data.  An  extensive  search  of  the  unclassified  literature  was  under¬ 
taken  to  recover  the  written  reports  containing  the  recorded  atmospheric  conditions  for 
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Table  5:  Bubble  Rise  Validation  Set  Sorted  by  Yield 


Operation 

Shot  Name 

Yield  (kt) 

HOB  (m) 

Hardtack  II 

Humboldt 

0.0078 

7.62 

Hardtack  n 

Catron 

0.021 

22.10 

Hardtack  II 

Vesta 

0.024 

0.00 

Hardtack  II 

DonaAna 

0.037 

137.16 

Hardtack  II 

Hidalgo 

0.077 

114.91 

Hardtack  II 

Quay 

0.079 

30.48 

Hardtack  II 

Eddy 

0.083 

152.40 

Hardtack  II 

RioAniba 

0.09 

22.10 

Hardtack  II 

Wrangell 

0.115 

457.20 

Plumbbob 

Franklin 

0.14 

91.44 

Piumbbob 

Wheeler 

0.197 

152.40 

Upshot-Knothole 

Ray 

0.2 

30.48 

Upshot-Knothole 

Ruth 

0.2 

92.87 

Sunbeam 

JohnnieBoy 

0.5 

-0.58 

Plumbbob 

Laplace 

1 

228.60 

Hardtack  n 

SantaFe 

1.3 

457.20 

Hardtack  II 

Lea 

1.4 

457.20 

Plumbbob 

John 

2 

6096.00 

Hardtack  n 

Mora 

2 

457.20 

Hardtack  II 

DeBaca 

2.2 

457.20 

Plumbbob 

FranklinPrime 

4.7 

228.60 

Hardtack  II 

Sanford 

4.9 

457.20 

Hardtack  n 

Socorro 

6 

441.96 

Plumbbob 

Morgan 

8 

152.40 

Plumbbob 

Owens 

9.7 

152.40 

Plumbbob 

Kepler 

10 

152.40 

Plumbbob 

Wilson 

10 

152.40 

Upshot-Knothole 

Dixie 

11 

1835.51 

Plumbbob 

Doppler 

11 

457.20 

Plumbbob 

Fizeau 

11 

152.40 

Plumbbob 

Galileo 

11 

152.40 

Plumbbob 

Boltzman 

12 

152.40 

Plumbbob 

Charleston 

12 

457.20 
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Table  5:  Bubble  Rise  Validation  Set  Sorted  by  Yield 


Operation 

Shot  Name 

Yield  (kt) 

HOB  (m) 

Plumbbob 

Newton 

12 

457.20 

Upshot-Knothole 

Grable 

15 

159.72 

Upshot-Knothole 

Annie 

16 

91.44 

Plumbbob 

Diablo 

17 

152.40 

Plumbbob 

Shasta 

17 

152.40 

Plumbbob 

Stokes 

19 

457.20 

Plumbbob 

Whitney 

19 

152.40 

Upshot-Knothole 

Badger 

23 

91.44 

Upshot-Knothole 

Nancy 

24 

91.44 

Upshot-Knothole 

Encore 

27 

738.53 

Upshot-Knothole 

Harry 

32 

91.44 

Plumbbob 

Priscilla 

37 

213.36 

Redwing 

Lacrosse 

40 

5.18 

Upshot-Knothole 

Simon 

43 

91.44 

Plumbbob 

Smoky 

44 

2)3.36 

Upshot-Knothole 

Climax 

6i 

406.60 

Plumbbob 

Hood 

74 

457.20 

Castle 

Koon 

110 

4.15 

Redwing 

Zuni 

3500 

2.74 

Redwing 

Tewa 

5000 

4.57 

Castle 

Bravo 

15000 

2.13 

these  shots  (26;27;40).  While  most  were  available  from  DAS  LAC  (8),  the  Plumbbob  series 
had  to  be  recovered  from  the  original  adiabatic  chart  entries.  Microfilm  recordings  of  these 
charts  are  held  at  Asheville,  NC,  and  copies  were  obtained  by  the  author  (50).  The  weather 
data  included  the  temperature,  pressure,  relative  humidity  (or  dew  point  temperature), 
wind  speed  and  direction  for  various  altitudes  from  the  surface  to  altitudes  above  the  stabi- 
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lized  cloud  top.  All  data  was  entered  into  the  computer  by  the  author  for  use  in  the  CRM 
code.  Every  effort  was  made  to  ensure  accuracy  in  the  data  entry. 

Eractional  Root  Mean.  SauareJERMS)  Error  definition.  The  figures  of  merit  used 
by  Norment  to  judge  the  performance  of  the  1977  version  of  the  CRM  were  the  FRMS 
deviation  and  the  fractional  mean  deviation  (FMD)  (46:22) 


FRMS  = 


(22) 
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(23) 


FMD 


zobs  ~  zcalc 

£ 

V.  Z0bs  ; 

N 

where 

2obs=°bserved  cloud  top  height  (m) 

2ca]c=calculated  cloud  top  height  (m) 

N  -  number  of  shots  in  data  set  (54) 

The  same  figures  of  merit  were  used  in  the  validation  of  changes  made  in  this  research. 
Relative  or  fractional  deviations  were  used  due  to  the  range  of  observed  cloud  heights 
used  in  these  studies. 
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The  purpose  of  this  chapter  is  to  explain  the  theory  of  particle  rise  in  DELFIC’s  CRM. 
Both  the  current  and  proposed  methods  of  particle  rise  are  dependent  on  the  results  of  the 
previous  chapter  (See  “Bubble  Rise  Dynamics”  on  page  11).  The  present  chapter  begins  in 
the  first  section  by  describing  the  levels  of  effort  for  panicle  rise  simulation.  The  second 
section  expands  on  the  method  of  using  vortex  theory  to  model  the  flow  of  particles.  The 
third  section  goes  into  the  details  of  how  this  method  is  implemented  in  DELFIC’s  CRM. 
The  final  section  reports  the  cases  that  were  run  and  the  method  of  comparison  to  hydro¬ 
code  data  and  the  current  particle  rise  of  the  CRM.  Actual  results  and  comparisons  will  be 
found  in  Chapter  IV. 

Levels  of  Effort 

There  are  three  main  ways  of  modeling  particle  rise  dynamics.  The  first  is  the  current 
method  employed  by  DELFIC’s  CRM  (one-dimensional).  In  addition,  there  are  two  pri¬ 
mary  ways  of  modeling  the  particle  flow  fields  in  two  dimensions.  The  first  is  to  use  a  2-D 
hydrodynamics  model.  The  second  is  to  use  an  analytic  solution  to  the  fluid  flow  for  a 
bubble  rising  in  a  still  fluid.  The  flow  field  method  is  a  computationally  cheaper  way  of 
obtaining  increased  fidelity  of  particle  rise  simulation  in  DELFIC  without  resorting  to 
computationally  expensive  hydrocodes. 

Current  DELFIC  Modeling.  This  is  a  simplified  description  of  a  very  complex  algo¬ 
rithm  used  to  determine  the  particle  distribution  in  the  stabilized  cloud.  The  1-D  particle 
rise  subroutine  in  the  CRM  takes  the  history  of  the  stored  bubble  rise  variables  and  traces 
the  upward  motion  of  an  initially  uniformly  distributed  cloud.  Then,  allowing  for  fallout 
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(or  particle  settling)  during  the  rise,  the  particles  are  followed  in  their  path  up  through  the 
atmosphere  to  stabilization.  Once  the  end  time  of  the  stored  variable  table  is  reached  (See 
Table  4  on  page  26),  the  particles  are  fixed  in  their  vertical  locations  and  the  stabilized 
cloud  is  thus  defined.  This  algorithm  then  stores  the  description  of  the  stabilized  cloud  so 
that  it  can  be  used  in  the  DTM  of  DELFIC  (47:27-30). 

Source  Initial  Location.  The  particle  cloud  is  represented  as  a  cylinder  with  a  top, 
bottom,  and  radius  that  match  the  top,  bottom,  and  radius  of  the  spheroidal  bubble  defined 
in  the  previous  chapter  (See  Figure  6  on  page  42).  The  mass  of  particles  in  the  cloud  is 
assumed  to  be  uniform  in  its  spatial  distribution  at  the  beginning  of  the  particle  rise  calcu¬ 
lations. 

Since  the  CRM  does  not  calculate  a  radial  distribution  of  mass,  the  cloud  is  only  sub¬ 
divided  into  a  set  of  vertical  layers.  The  number  of  layers,  KDI,  is  input  by  the  user  or  cho¬ 
sen  by  the  code  to  be  KDI  =  15  +  In  ( W  in  kilotons).  These  uniformly  loaded  layers  are 
represented  as  cylindrical  disks  (or  parcels).  This  vertical  stack  of  disks  in  the  initial  cloud 
is  repeated  for  each  particle  size  bin  in  the  particle  size  distribution.  The  particle  size  bins 
are  represented  by  a  geometric  mean  group  diameter,  bounded  by  upper  and  lower  bound¬ 
ary  particle  diameters. 

Particle  Size  Distribution.  DELFIC  contains  three  different  particle  size  distribu¬ 
tion  options:  log  normal,  power  law,  and  tabular  input.  The  default  distribution  chosen  by 
the  code  is  log  normal.  This  log  normal  distribution  has  a  mean  diameter  of  130  microns 
and  a  geometric  standard  deviation  of  4.  This  distribution  was  determin° 4 '  om  downwind 
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Figure  6:  Cross  Section  of  Current  CRM  Initial  Particle  Cylinder  With  Disk 
Structure  Based  on  Bubble  Vertical  and  Radial  Dimensions  (dashed)  in  Meters 

samples  of  fallout  and  is  well  accepted  by  the  fallout  community.  The  particle  sizes  for  a 
50  group  distribution  are  shown  in  Table  6  on  page  43. 

The  hydrocode  TASS  on  the  other  hand  has  two  different  particle  size  distributions: 
one  for  ejecta  mass  and  one  for  scoured  mass.  Both  distributions  are  determined  using 
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Table  6:  Default  Particle  Size  Distribution  for  DELFIC  (50  Groups) 


up 

Mean  Diameter  (m) 

Lower  Boundary  (m) 

Upper  Boundary  (m) 

1 

0.32686E-02 

0.22399B-02 

0.47695E-02 

2 

0.18155E-02 

0.14714E-02 

0.22399E-02 

3 

0.12845E-02 

0.11214E-02 

0.14714E-02 

4 

0.10108E-02 

0.91113E-03 

0.11214E-02 

5 

0.83632E-03 

0.76766E-03 

0.91113E-03 

6 

0.71296E-03 

0.662 17E-03 

0.76766E-03 

7 

0.62008E-03 

0.58067E-03 

0.66217E-03 

8 

0.54709E-03 

0.51546E-03 

0.58067B-03 

9 

0.48794E-03 

0.46189E-03 

0.5 1546E-03 

10 

0.43885E-03 

0.41696E-03 

0.46189E-03 

11 

0.39736E-03 

0.37867E-03 

0.41696E-03 

12 

0.36175E-03 

0.34559E-03 

0.37867E-03 

13 

0.33082E-03 

0.31669E-03 

0.34559E-03 

14 

0.30367E-03 

0.29119E-03 

0.31669E-03 

15 

0.27962E-03 

0.26852E-03 

0.291 19E-03 

16 

0.258 16E-03 

0.2482  IE-03 

0.26852E-03 

17 

0.23888E-03 

0.22991E-03 

0.24821E-03 

18 

0.22146E-03 

0.21332E-03 

0.22991E-03 

19 

0.20563E-03 

0. 19822E-03 

0.21332E-03 

20 

0.191 19E-03 

0.18440E-03 

0.19822E-03 

21 

0.17795E-03 

0.17172E-03 

0.18440E-03 

22 

0.16577E-03 

0.16003E-03 

0.17172E-03 

23 

0.15453E-03 

0.14922E-03 

0.16003E-03 

24 

0.14412E-03 

0.13919E-03 

0.14922E-03 

25 

0.13445E-03 

0.12986E-03 

0.13919E-03 

26 

0.12544E-03 

0.12116E-03 

0.12986E-03 

27 

0.11702E-03 

0.11302E-03 

0.12116E-03 

28 

0.10914E-03 

0.10539E-03 

0.11302E-03 

29 

0.10174E-03 

0.982 12E-04 

0.10539E-03 

30 

0.94774E-04 

0.91456E-04 

0.98212E-04 

31 

0.882 11E-04 

0.85082E-04 

0.91456E-04 

32 

0.820 15E-04 

0.79058E-04 

0.85082E-04 

33 

0.76153E-04 

0.73355E-04 

0.79058E-04 
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Table  6:  Default  Particle  Size  Distribution  for  DRLFIC  (50  Groups) 


Group 

Mean  Diameter  (m) 

Lower  Boundary  (m) 

Upper  Boundary  (m) 

34 

0.70599E-04 

0.67946E-04 

0.73355E-04 

35 

0.65326E-04 

0.62807E-04 

0.67946E-04 

36 

0.603 13E-04 

0.579 17E-04 

0.62807E-04 

37 

0.55537E-04 

0.53254E-04 

0.579 17E-04 

38 

0.50979E-04 

0.48800E-04 

0.53254E-04 

39 

0.46620E-04 

0.44537E-04 

0.48800E-04 

40 

0.42443E-04 

0.40447E-04 

0.44537E-04 

41 

0.38430E-04 

0.365 13E-04 

0.40447E-04 

42 

0.34564E-04 

0.327 18E-04 

0.365 13E-04 

43 

0.30826E-04 

0.29044E-04 

0.327 18E-04 

44 

0.27 198E-04 

0.25469E-04 

0.29044E-04 

45 

0.23655E-04 

0.21969E-04 

0.25469E-04 

46 

0.20165E-04 

0.18510E-04 

0.21969E-04 

47 

0.16685E-04 

0.15039E-04 

0.18510E-04 

48 

0.13129E-04 

0.11461E-04 

0.15039E-04 

49 

0.92895E-05 

0.75292E-05 

0.11461E-04 

50 

0.51597E-05 

0.35360E-05 

0.75292E-05 

power  law  relationships  and  range  in  size  from  1  micron  to  2  centimeters  for  scoured  mass 
and  50  microns  to  1  meter  for  ejecta  mass  (9:3).  Since  the  hydrocodes  use  dynamic  load¬ 
ing  for  assigning  mass  from  these  two  phenomena,  it  cannot  be  determined  a  priori  what 
the  particle  size  distribution  of  the  total  cloud  mass  will  be  at  stabilization.  For  this  reason, 
and  the  fact  that  the  particle  size  distribution  of  DELF1C  is  based  on  fallout  studies,  the 
default  distribution  of  DELFIC  is  used  in  the  current  research. 

One-Dimensional  Particle  Rise.  An  individual  parcel  (a  given  particle  size  at  a 
given  altitude)  is  traced  through  the  bubble  rise  variable  table  (See  Table  4  - page  26)  to 
determine  its  rise  history.  For  particle  rise  simulations,  the  top  and  base  of  each  parcel  are 
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followed  separately  with  the  upper  particle  size  boundary  corresponding  to  the  base  of  the 
parcel  and  the  lower  particle  size  boundary  used  for  the  parcel  top.  The  upward  velocity  of 
the  top  or  base  is  determined  by  its  location  relative  to  the  bubble’s  top  and  base  locations 
for  the  particular  time  in  its  history. 

For  a  given  time,  the  upward  velocity  for  a  parcel  wafer  (top  or  base)  is  found  from  the 
top  and  base  rise  rates.  If  the  wafer  is  within  the  main  cloud,  the  velocity  is  found  by  lin¬ 
early  interpolating  with  altitude  between  the  top  and  base  velocity  values.  A  wafer  exactly 
at  the  cloud  center  altitude  will  have  an  upward  velocity  equal  to  the  average  of  the  top 
and  base  velocities.  If  the  wafer  is  below  the  cloud  base,  the  velocity  is  calculated  by  lin¬ 
early  interpolating  between  the  base  rise  rate  and  a  zero  rise  rate  for  the  surface  altitude. 

Gravity  Settling.  The  only  way  a  wafer  that  starts  within  the  main  cloud  can  get 
below  the  cloud  base  is  to  have  a  downward  component  to  its  velocity  vector.  This  compo¬ 
nent  comes  from  gravity  settling.  Larger  particles  will  have  a  higher  terminal  velocity  due 
to  gravity.  This  implies  that  the  base  wafers  represented  by  the  upper  boundary  size  will 
fall  faster  than  the  top  wafers  which  use  the  lower  boundary  particle  size.  This,  combined 
with  the  fact  that  the  cloud  does  not  collapse  but  rather  expands,  means  that  individual 
parcels  are  stretched  vertically. 

If  a  wafer  never  leaves  the  main  cloud,  its  final  radius  will  be  that  of  the  main  cloud  at 
vertical  stabilization.  If  wafers  fall  below  the  base  of  the  main  cloud,  their  radius  is  set  to 
the  radius  of  the  cloud  at  the  time  it  departs  the  base  of  the  main  cloud.  Once  a  parcel  has 
reached  the  end  of  the  bubble  rise  history,  it  stops  its  vertical  rise.  If  a  parcel  base  and  top 
have  different  radii,  the  parcel  may  be  subdivided  in  the  vertical  dimension.  Each  subdi¬ 
vided  parcel  is  given  a  radius  depending  on  its  location  in  the  tapered  parent  parcel.  This 
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method  of  subdividing  parcels  which  have  at  least  partially  fallen  out  of  the  main  cloud  is 
how  a  fallout  stem  is  defined  in  the  DELFIC  stabilized  cloud. 

Shot  Time  Wind  Shifting.  The  final  subroutine  in  the  CRM  adjusts  the  stabilized 
parcels  in  the  horizontal  direction  due  to  the  ambient  winds.  This  is  a  way  of  modeling  the 
asymmetry  of  stabilized  clouds  above  ground  zero  that  were  observed  in  the  atmospheric 
nuclear  tests.  Each  parcel’s  center  is  adjusted  according  to  the  wind  velocity  and  the 
amount  of  time  it  spends  in  each  wind  altitude  layer. 

Original  CRTTM  Included  a  2-D  Option.  As  a  note,  it  should  be  pointed  out  that 
the  1967  version  of  the  Cloud  Rise^Transport  Interface  Module  (CRTTM)  had  the  option  of 
accepting  input  for  a  two-dimensional  description  of  the  particle  cloud  (22:3-1,3-2).  This 
capability  was  in  anticipation  of  someday  linking  its  transport  module  to  a  hydrocode 
solution  of  the  rising  cloud.  Unfortunately  this  capability  was  removed  from  the  1-D  parti¬ 
cle  rise  and  wind  shifting  subroutines  which  describe  the  stabilized  particle  cloud. 

Modeling  with  Hydrocodes.  A  second  method  of  computing  the  particle’s  rise 
dynamics  is  with  hydrocodes.  Hydrocodes  follow  the  variables  in  a  multicell  framework. 
This  allows  the  variables  to  change  in  multiple  dimensions.  As  the  result  of  solving  a  set 
of  partial  differential  equations,  a  two-dimensional  representation  of  all  the  variable  solu¬ 
tions  are  determined.  Therefore  the  velocity  of  the  flow  is  known  with  both  vertical  and 
radial  components. 

With  hydrocodes,  the  stabilized  dust  cloud  is  determined  by  the  use  of  trace  particles. 
Each  trace  particle  represents  a  certain  amount  of  mass  in  the  cloud.  Trace  particles  are 
added  when  modeling  the  addition  of  mass  from  ejecta  or  sweepup.  The  sizes  and  weights 
of  these  trace  particles  are  determined  by  the  particle  size  distributions  used  in  these  codes 
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and  the  amount  of  soil  being  injected  into  the  cloud.  Once  in  the  multicell  framework,  the 
particles  are  moved  by  the  forces  determined  by  the  local  state  of  the  cloud. 

Modeling  with  Vortex  Theory.  An  alternative  to  solving  all  of  the  variables  in  two 
dimensions  is  to  use  the  theory  of  bubbles  rising  in  a  still  fluid.  The  conservation  equa¬ 
tions  as  described  in  the  previous  chapter  can  be  solved  and  the  bubble’s  rise  history 
stored.  By  knowing  just  the  rise  velocity  and  dimensions  of  the  bubble,  the  two-dimen¬ 
sional  flow  in  and  around  the  main  cloud  can  be  analytically  determined.  The  vortex  flow 
field  models  the  circulation  within  the  rising  cloud.  Some  assumptions  concerning  the 
fluid  flow  are  needed.  The  current  1-D  CRM  approach  contains  no  circulation! 

VORDUM.  One  such  code  which  used  this  approach  is  the  Vortex  Dust  Model 
(VORDUM).  VORDUM  uses  empirical  fits  to  the  velocity  and  dimensions  of  the  bubble 
as  a  function  of  time  instead  of  solving  a  set  of  conservation  equations  (7;64;65).  Its  fits, 
therefore,  do  not  take  into  account  any  variation  due  to  atmospheric  conditions  in  comput¬ 
ing  the  bubble  rise  solution.  The  empiricisms  are  based  solely  on  yield  and  height  of  burst. 

Radioactive  Dust.  DELFIC’s  CRM  is  concerned  solely  with  the  radiation  hazard 
of  the  fission  fragments  as  they  appear  in  the  engulfed  soil  and  weapon  debris  mixture. 
VORDUM  is  a  systems-level  code  for  producing  the  total  dust  environment  (radioactive 
and  nonradioactive)  of  a  nuclear  explosion  for  fratricide  and  radar  communication  studies. 
It  makes  no  attempt  to  determine  what  dust  is  radioactive,  nor  does  it  concern  itself  with 
dust  particles  once  they  are  grounded.  This  is  in  contrast  to  DELFIC  which  ultimately  cal¬ 
culates  the  location  and  intensity  of  only  the  grounded  radioactivity. 

Circulation.  A  limitation  to  the  vortex  flow  theory,  is  that  the  flow  is  directly 
related  to  the  rise  velocity.  This  means  that  when  the  bubble  stops  rising  vertically,  the 
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model  will  revert  back  to  a  one-dimensional  gravity  settling  model  (equivalent  to  the  cur¬ 
rent  CRM).  If  only  from  flow  inertia,  one  would  expect  some  residual  circulation  after  the 
rise  has  stopped.  To  model  the  residual  circulation  after  the  rise  has  stopped,  a  full  2-D 
hydrodynamics  solution  could  be  used,  but  not  without  the  added  expense  of  running  such 
a  code. 

Theory  Used  in  Vortex  Particle  Rise  Module 

The  motivation  for  using  an  analytic  flow  field  solution  is  that  it  matches  the  observed 
patterns  of  rising  nuclear  clouds.  The  toroidal  shape  of  the  flow  pattern  (See  Figure  1  on 
page  3)  is  reproduced  with  such  a  solution,  unlike  the  current  DELFIC  modeling.  The 
cloud  properties  needed  for  the  solution  are  available  from  the  bubble  rise  solution. 
Finally,  the  velocity  distribution  predicted  by  this  type  of  solution  has  been  shown  to 
match  laboratory  experiments  (30;  60). 

To  understand  vortex  flow  fields,  some  theory  from  classical  fluid  dynamics  needs  to 
be  reviewed.  The  topics  covered  in  this  section  begin  with  a  discussion  of  vortices  in  gen¬ 
eral  and  Hill’s  spherical  vortex  in  particular.  From  this,  the  definition  of  a  streamline  is 
presented  with  application  to  the  rising  bubble.  Velocity  vectors  can  then  be  derived  from 
the  knowledge  of  the  streamlines,  and  used  for  particle  motion.  The  application  of  an 
image  vortex  is  presented  which  will  help  account  for  the  boundary  with  the  ground. 
Some  of  the  more  recent  theory  covering  spheroidal  vortices  is  then  discussed. 

Vortices.  When  a  buoyant  bubble  rises  in  a  denser  atmosphere  a  vortex  circulation  is 
induced.  Analytic  solutions  to  this  circulation  are  found  within  classical  fluid  dynamics 
texts  (3:526;  19;  29:245;  37:578-579;  57).  There  are  three  main  categories  for  these  solu- 
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dons:  ling  vortices,  spherical  vortices,  and  spheroidal  vortices.  Ring  vortices  differ  from 
the  other  two  types  in  that  the  vorticity  is  concentrated  in  a  smaller  core  region  of  the  flow 
(56:7). 

When  modeling  nuclear  clouds  with  analytic  solutions,  certain  assumptions  are 
required.  The  assumptions  which  allow  one  to  use  an  analytic  solution  to  the  flow  field  are 
as  follows.  It  is  assumed  that  the  flow  in  and  around  the  rising  bubble  is  steady.  Also,  it  is 
assumed  that  the  fluid  is  inviscid  and  incompressible.  The  irrotational  flow  outside  the 
bubble  is  assumed  to  be  the  potential  flow  in  a  still  fluid  about  a  moving  sphere.  The  nor¬ 
mal  component  of  the  velocity  must  vanish  at  the  sphere  boundary  while  the  tangential 
velocity  must  remain  continuous  (56:12). 

Hill’s  Spherical  Vortex  .  Hill’s  spherical  vortex  is  described  by  Batchelor  as  “one 
extreme  member  of  a  family  of  vortex  rings”  (3:526).  The  extremity  is  that  its  vorticity 
occupies  the  entire  sphere,  rather  than  only  a  limited  region. 

The  solution  for  Hill’s  spherical  vortex  is  obtained  by  solving  for  the  flow  past  a  sta¬ 
tionary  sphere.  The  solution  for  a  sphere  moving  in  a  still  fluid  is  then  obtained  by  super¬ 
imposing  on  the  entire  system  an  equal  but  opposite  potential  flow.  Hill  showed  that  a 
solution  to  the  equations  below  predicts  the  flow  in  and  about  a  fluid  moving  past  a  sta¬ 
tionary  sphere  (56:13). 


1 

r2(sin0)2 


=  0 


for  r  <,  a 


(24) 
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(25) 


1  I’d  f  1  dy)  d  f  1  dy^l 

r2  ( sin0) 2  sin0dr  /  30Vr2sin05F' J 


=  -4  (y) 


for  r>a 


where 

r  »  radius  coordinate  for  spherical  coordinate  system  (m) 

0  =  angle  from  axis  of  motion  in  spherical  coordinate  system 
V  =  streamfunction  of  the  flow  to  be  solved  for 
a  *  radius  of  the  sphere  (m) 

A(V)=  arbitrary  function  of  y 

Streamlines  and  Streamfunctions.  A  streamline  is  a  curve  everywhere  parallel  to  the 
direction  of  flow.  Streamlines  represent  lines  of  constant  streamfunction,  y  (53:23-24). 
The  above  equations  are  solved  by  requiring  the  streamfunction  to  vanish,  and  the  tangen¬ 
tial  velocity  to  be  continuous,  at  the  sphere  surface,  i.e., 

V  =  0  at  r  -  a  (26) 

rslnfldr'  **  cont*nuous  at  r  ~  a  (27) 

The  final  boundary  condition  is  that  the  flow  approaches  the  potential  flow  about  a  solid 
sphere  as  it  approaches  the  free  field  stream  (56: 13). 
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(28) 


\y  — »  --t/r2 ( sin0) 2  as  r— »» 


where 

U  =  velocity  of  the  moving  body  of  fluid  (m/s) 

The  above  equations  and  boundary  conditions  are  solved  to  obtain  the  solution  for  a 

1  2 

fluid  flowing  past  a  stationary  sphere.  The  free  field  stream  function  ( -  Ur 2  ( sin0)  )  is 

superimposed  on  the  previous  answer  to  get  the  solution  for  a  moving  sphere  in  a  still 
fluid.  The  streamfunctions  for  the  flow  inside  and  outside  Hill’s  spherical  vortex  moving 
in  a  still  fluid  are  then  (56: 175) 


vy  s  Ur2  (  sin0) j*  -  inside  sphere  (29) 

2 

outside  sphere  (30) 

The  streamlines  for  this  solution  are  shown  in  Figure  7  on  page  52. 

Velocity  Vectors.  The  velocity  components  for  both  inside  and  outside  the  spherical 
vortex  can  be  found  by  differentiating  the  appropriate  streamfunction.  The  radial  and  tan¬ 
gential  velocities  are  defined  by  (56:35-36) 
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Vertical  Direction 


Figure  7:  Streamlines  for  Hill’s  Spherical  Vortex  (cross  sectional  view  of  unit  sphere) 


Ur  = 


1  dy 
r2sin  0^9 


(31) 


u  = _ 

0  rsinQdr 


(32) 
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This  yields  a  radial  velocity  of 


Ur  =  t/co$e(l-^) 
ar  =  -|t/c0Se(  i  -  £|) 

and  a  tangential  velocity  of 

Ue  =  -£/sin0^1  + 

t/e  =  |t/sine(  1  - 


outside  sphere 


inside  sphere 


outside  sphere 


inside  sphere 


(33) 


(34) 


(35) 


(36) 


These  equations  are  then  transformed  into  cylindrical  coordinates  for  implementation  in 
the  code  to  give  the  vertical  and  radial  components  of  the  particle’s  velocity.  A  vector  field 
of  the  velocity  shows  the  magnitude  of  the  flow  in  different  locations,  whereas  the  stream¬ 
line  figure  gave  only  a  qualitative  view  of  the  direction  of  flow  (See  Figure  8  on  page  54). 

From  this  figure,  one  can  see  the  difference  compared  to  the  1-D  model  of  the  current 
CRM.  Not  only  is  there  the  addition  of  a  radial  component,  but  there  is  variation  of  the 
vertical  component  based  on  location  in  the  sphere.  The  vertical  component  varies  from 
2.5U  at  the  center  to  -0.5U  at  the  edge  of  the  sphere,  along  the  vortex  midplane. 
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Figure  8:  Vector  Field  for  Hill’s  Spherical  Vortex  (unit  sphere  with  unit  velocity) 


Image  Vortex.  A  mirror  image  of  the  spherical  vortex  is  created  at  an  equal  distance 
on  the  other  side  of  the  ground  plane.  This  new  vortex  has  the  same  dimensions  and  speed 
as  the  real  vortex  but  moves  in  the  opposite  direction.  This  has.  the  effect  of  canceling  out 


the  vertical  components  and  doubling  the  radial  component  at  the  ground  surface.  This  not 
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Figure  9:  Direction  and  Location  of  Image  Vortex  Sphere  (35:42) 


i 
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only  satisfies  the  condition  that  flow  cannot  enter  the  ground  plane,  but  also  helps  simulate 
the  intense  surface  winds  which  are  known  to  occur  for  near  surface  nuclear  bursts 
(56:60). 

After  the  image  vortex  is  included,  there  are  now  three  possible  locations  for  a  parti¬ 
cle:  inside  both  the  real  and  imaginary  vortices,  inside  the  real  vortex  and  outside  the 
imaginary  vortex,  and  outside  both  the  real  and  imaginary  vortices.  After  conversion  to 
cylindrical  coordinates  and  inclusion  of  the  image’s  effect,  the  vertical  and  radial  compo¬ 
nents  of  velocity  are  found  to  be  the  following. 

Inside  both  the  real  and  imaginary  vortices: 
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where 

Uz  =  the  vertical  velocity  (m/s) 

UR  *  the  radial  velocity  (m/s) 

R  =  the  radial  distance  from  the  axis  of  symmetry  (m) 

Zj  =  the  vertical  distance  from  the  midplane  of  the  real  vortex  (m) 

Z2  -  the  vertical  distance  from  the  midplane  of  the  imaginary  vortex  (m) 
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Inside  the  real  vortex  and  outside  the  imaginary  vortex: 
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Outside  both  the  real  and  imaginary  vortices: 
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The  effect  of  the  image  vortex  on  the  velocity  field  for  a  sphere  intersecting  the  ground  is 
shown  in  Figure  10  on  page  58.  Note  the  vectors  at  ground  level  have  no  vertical  compo¬ 
nent. 

Spheroidal  Vortices.  A  solution  similar  to  Hill’s  spherical  vortex  was  found  for  sphe¬ 
roids  by  O’Brien  (49)  This  solution  was  investigated  by  this  and  other  authors,  but  not 
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Figure  10:  Velocity  Field  for  Sphere  Intersecting  the  Ground  (Image  Vortex 

Included) 


implemented  in  the  current  research.  It  was  implemented  in  the  VORDUM  dust  code  by 
Shannon  and  he  reports  that  “results  from  the  spheroidal  vortex  VORDUM  model  do  not 
differ  greatly  from  the  spherical  vortex  VORDUM  model  results.,. and  applications  of  the 


58 


spheroidal  vortex  flow  field  to  define  dust  environments  for  use  in  system  studies  may  not 
be  warranted”  (56: 174). 

Perturbation  Theory.  A  final  note  concerning  the  results  of  this  author’s  literature 
search  is  in  order.  Pozrikidis  mathematically  showed  that  an  initially  spherical  vortex 
which  becomes  transformed  (due  to  axisymmetric  perturbations)  to  an  oblate  shape, 
entrains  the  outside  irrotationa!  fluid.  This  fluid  enters  the  rotational  fluid  from  the  rear 
stagnation  point  and  travels  along  the  vortex  boundary.  He  states  that  this  action  reduces 
the  oblate  vortex  to  a  nearly  steady  vortex  ring  (51). 

His  mathematical  description  is  qualitatively  similar  to  the  shape  history  of  real 
nuclear  clouds.  Such  clouds  start  as  near  perfect  spherical  fireballs  and  are  deformed  dur¬ 
ing  their  intense  motion  upwards  into  oblate  spheroids.  Some  researchers  have  even  mod¬ 
eled  nuclear  clouds  as  vortex  rings  (11). 

Vortex  Implementation  in  DELFIC’s  Current  CRM 

On  the  basis  of  the  previous  theory,  this  section  describes  how  the  author  has  imple¬ 
mented  vortex  flow  into  the  CRM  of  DELFIC  to  produce  a  two-dimensional  description  of 
the  stabilized  cloud.  The  major  subsections  discuss  the  modeling  choices  made  by  the 
author  to  complement  the  theory  of  the  previous  section.  The  first  subsection  shows  how 
the  bubble  rise  history  is  used  in  the  description  of  the  vortex  region.  The  second  subsec¬ 
tion  reviews  the  modeling  of  the  initial  soil  distribution  for  the  radioactive  dust.  This  is 
followed  by  a  brief  review  of  particle  settling  and  the  force  balance  needed  to  determine 
the  particle  acceleration  equations.  The  last  major  phenomenon  modeled  is  the  correction 
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for  the  expansion  of  the  bubble.  The  final  c"bsection  describes  the  solution  method  used  to 
calculate  the  particle  histories. 

Input  from  Bubble  Rise  History.  As  in  the  current  CRM,  particle  rise  is  dependent  on 
the  results  of  the  bubble  rise  equations.  The  history  of  the  cloud’s  velocity  and  dimensions 
are  needed  to  describe  fully  the  flow  field  using  the  theory  described  above.  In  addition, 
the  temperature  of  the  cloud  and  its  density  are  needed  to  determine  the  viscosity  and  the 
drag  force  on  the  particles.  A  major  question  that  still  remains  is  how  one  relates  the  visi¬ 
ble  oblate  bubble  to  the  domain  of  the  vortex  sphere. 

Visible  Cloud  Spheroid  vs.  Vortex  Sphere.  In  matching  the  oblate  spheroid  to  a 
domain  more  conducive  to  flow  field  simulation,  some  modeling  choices  had  to  be  made. 
In  the  present  research  the  top  of  the  bubble  was  matched  to  the  top  of  the  vortex  sphere 
region.  The  radius  of  the  sphere  was  then  matched  to  the  radius  of  the  bubble.  This  pro¬ 
duced  a  vortex  region  which  encompassed  most  of  the  visible  cloud  and  included  a  region 
below  the  cloud  (See  Figure  11  on  page  61). 

The  region  below  the  cloud  but  still  within  the  vortex  sphere  helps  simulate  the  updraft 
in  the  stem.  Also,  since  the  stagnation  point  for  the  center  of  circulation  of  a  rising  sphere 
is  above  the  midplane  of  the  expanding  sphere,  the  particles  will  concentrate  within  the 
visible  oblate  sphere  which  is  also  centered  above  the  midplane  (60:201). 

Initial  Particle  Locations.  The  initial  location  of  the  soil  in  the  cloud  must  now  be 
determined  and  a  set  of  trace  particles  described.  Each  of  these  trace  particles  must  be 
assigned  a  fraction  of  the  soil  mass  based  on  their  location  and  an  assumption  of  the  soil 
loading  in  the  cloud  before  the  cloud  begins  its  rise.  The  effects  of  sweepup  and  time 
dependent  loading  will  also  be  discussed. 
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Main  Source  Cylinder.  The  initial  location  of  the  soil  particles  is  assumed  to  be  a 
cylinder  with  the  initial  cloud  top  and  radius  (See  Figure  12  on  page  62).  This  is  similar  to 
the  current  assumption  in  DELFIC’s  CRM  with  one  difference.  The  difference  in  the  new 
model  is  that  the  cylinder  is  assumed  to  be  touching  the  ground  instead  of  hovering  in 
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space.  This  accounts  for  the  pedestal  of  very  early  time  dust  clouds  and  is  similar  to  the 
way  VORDUM  originally  modeled  its  ejecta  source  (7:204). 


Figure  12:  Vortex  Particle  Source  Cylinder  and  Initial  Bubble  (both  dashed)  and 
VORDUM’s  2-Cylinder  Sources  in  Relation  to  Initial  Fireball  (distances  in  meters) 


Sweepup  Source  Cylinder.  The  original  VORDUM  also  included  a  second  source 
cylinder,  a  shorter  but  wider  sweepup  cylinder.  This  second  cylinder  modeled  the  dust 
which  is  not  directly  engulfed  by  the  fireball  but  is  eventually  swept  into  the  fireball  or  its 
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stem  due  to  the  flow-induced  surface  winds.  Although  VORDUM’s  dust  loading  was 
immediate,  as  is  DELFIC’s,  the  initial  time  of  particle  tracking  was  later  than  in  DELFIC 
(nominally  30  seconds).  The  question  remains  which  of  this  dust  is  part  of  the  radioactive 
dust  that  eventually  contributes  to  the  fallout  dose  on  the  ground. 

McGahan  showed  that  die  debris  in  the  ejecta  (VORDUM’s  first  cylinder)  has  a  vastly 
different  final  location  in  the  stabilized  cloud  when  compared  to  the  stabilized  sweepup 
dust.  One  may  conjecture  that  the  nonradioactive  sweepup  dust  may  attach  itself  to  the 
radioactive  ejecta  debris  if  they  are  in  the  same  location  of  the  cloud.  This  area  of  overlap 
between  the  two  sources  appears  to  be  small  from  initial  investigations  (36:77).  One 
source  cylinder  is  used  in  this  research. 

The  main  source  cylinder  dimensions  chosen  in  this  model  virtually  engulf  the  VOR- 
DUM  two  cylinder  source.  The  vortex  sphere  domain  in  DELFIC  also  starts  earlier  in  the 
rise,  and  therefore  at  a  lower  altitude.  The  effect  is  that  the  majority  of  the  particles  start 
above  the  vortex  midplane  and  with  a  greater  initial  radial  distribution.  It  was  found  that 
the  particles  which  start  lower  in  the  domain,  or  closer  to  the  axis,  take  a  longer  time  to 
enter  the  toroidal  core,  if  at  all.  The  particles  which  start  higher  in  the  domain,  or  away 
from  the  axis,  enter  the  toroidal  core  faster. 

Time  Dependent  Loading.  As  in  the  original  CRM,  the  new  model  for  particle  rise 
does  not  include  a  time  dependent  loading  of  the  cloud.  Instead,  the  entire  soil  mass  is 
assumed  to  be  in  the  cloud  at  the  initial  time.  The  soil  that  does  enter  at  later  times  (not 
modeled  here)  may  not  make  it  into  the  vortex  sphere  and  if  so  may  not  end  up  being  col¬ 
located  with  the  radioactive  debris  from  the  weapon. 
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Also  all  panicles  are  assumed  to  be  solid  from  the  start  with  essentially  their  final 
radius,  and  are  tracked  accordingly.  This  is  a  minor  approximation  since  all  clouds  simu¬ 
lated  in  this  study  cooled  to  below  soil  solidification  after  the  first  few  seconds  of  their  rise 
simulations  (</i  +  5  sec). 

Number  of  Trace  Particles.  In  determining  the  number  of  trace  particles,  one  must 
choose  a  number  small  enough  for  reasonable  run  times  but  large  enough  for  reasonable 
approximation  of  the  radioactive  dust  densities.  In  the  current  CRM,  the  cloud  is  divided 
into  a  one-dimensional  stack  of  cylindrical  parcel  To  make  the  new  option  for  particle 
rise  easy  to  implement,  this  same  number  of  vertical  grid  points  was  chosen.  In  addition, 
the  new  routine  has  KD1  horizontal  grid  points  for  starting  trace  particles.  This  gives  KDI 2 
number  of  trace  particles  for  each  size  class  run.  For  a  typical  run  with  50  size  classes,  this 
gives  approximately  20,000  (=20x20x50)  trace  particles  to  describe  the  soil  loading  of  the 
cloud  at  any  time.  This  number  met  both  of  the  aforementioned  criteria  and  had  the  ease  of 
user  input  already  set  up  in  the  current  CRM. 

Particle  Weighting.  Once  the  number  of  trace  particles  is  chosen,  their  starting 
locations  within  a  2-D  vertical  cross  section  of  the  initial  source  cylinder  are  determined. 
Each  particle  size  group  contains  an  identical  uniform  grid  of  trace  particles  from  the  axis 
to  the  out;r  cylinder  radius.  A  fraction  of  the  cylinder’s  volume  is  determined  by  the  3-D 
annulus  the  trace  particle  represents.  The  particle’s  weight  is  then  determined  by  the  num¬ 
ber  of  particle  sizes  and  the  fraction  of  the  cylinder’s  volume  its  starting  annulus  repre¬ 
sents.  This  amount  of  mass  weighting  stays  with  the  particle  and  is  output  with  the  history 
of  the  particle’s  rise  for  loading  information. 
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can  be  found  by  equating  the 


force  on  the  particle  to  the  sum  of  the  forces  of  gravity  and  drag  (35:48). 


f,’f,*fd 


(43) 


(VZ-Vz)[(Uz-V2)2+  (t/,-^)2]2  (44) 


where 

Fp  =  force  on  the  particle  (kg  m/s2) 
Fg  =  gravity  force  (kg  m/s2) 

Fd  =  drag  force  (kg  m/s2) 

Vg  =  particle  vertical  velocity  (m/s) 
VR  =  particle  radial  velocity  (m/s) 
rp  =  particle  radius  (m) 
d  -  particle  diameter  (m) 
pm  =  particle  (mass)  density  (kg/m3) 
pa  =  air  (or  cloud)  density  (kg/m3) 
Cp  =  coefficient  of  drag 
which  becomes 
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(45) 


v7  - 


l-;r7(fr  vz>  Oz -  vz) 2  +  (Uf  -  VV  *] 


4Pm 


1 

^ =  -IrT  (V«-£/«)[(c,z-vy2+ w 

An  alteraarive  to  calculating  the  above  accelerations  is  to  assume  that  the  particle  is  in 
equilibrium  with  the  flow  field.  The  particles  are  then  given  the  flow  velocity  plus  a  cor¬ 
rection  for  gravity.  This  correction  is  called  the  slip  correction.  The  slip  correction 
assumes  the  particle  immediately  achieves  its  terminal  velocity  of  fall.  The  method  of  cal¬ 
culating  this  terminal  velocity  is  the  same  as  in  the  current  CRM  (47:24-25).  The  particle’s 
location  is  determined  to  be  in  the  cloud  or  out  of  the  cloud,  and  its  corresponding  envi¬ 
ronment  temperature  and  density  are  found.  With  this  information  and  the  particle’s  diam¬ 
eter,  its  viscosity  and  drag  coefficient  can  be  calculated. 

Lagrangian  Transport  vs.  Slip  Correction.  The  hydrocode  TASS  uses  a 
Lagrangian  transport  routine  to  follow  particles  greater  than  1000  microns  (1  mm).  This  is 
in  contrast  to  assigning  the  particle  the  flow  velocity  minus  the  terminal  velocity  in  the 
vertical  direction.  This  assumes  that  particles  less  than  a  millimeter  are  in  equilibrium 
with  the  flow,  while  those  above  may  not  be.  This  criteria  was  maintained  in  the  current 
research. 

The  cutoff  value  use  d  in  TASS  (1  mm)  was  investigated  in  the  vortex  flow  field  code 
developed.  It  was  discovered  that  not  only  do  one  millimeter  particles  behave  as  if  they 
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are  in  equilibrium  with  the  flow,  but  so  do  particles  as  large  as  2.5  mm.  Particles  larger 
than  this  may  not  be  in  equilibrium  with  the  flow  but  their  fraction  of  the  particle  size  dis¬ 
tribution  for  DELFIC  is  negligible  (<2%).  Therefore  using  the  assumption  of  being  in 
equilibrium  with  the  flow  and  substituting  a  slip  velocity  for  Lagrangian  transport  will  not 
significantly  affect  the  result  of  DELFIC ’s  particle  rise.  The  savings  in  computer  time 
when  using  a  slip  correction  is  substantial,  and  therefore  the  switch  for  Lagrangian  trans¬ 
port  in  the  new  particle  rise  model  of  DELFIC  may  be  removed  if  desired. 

TASS  vs.  DELFIC.  The  DELFIC  settling  subroutine  is  based  on  the  work  of 
Davies  and  Beard  (47:25)  The  particle  settling  equations  used  in  TASS  are  based  on  a 
modified  Stokes  drag  term  (2: 19).  In  TASS,  the  drag  coefficient  for  a  perfect  sphere  is 
multiplied  by  1.5  for  a  rough  sphere  approximation.  With  this  method  of  calculating  the 
drag,  TASS  will  predict  a  slightly  larger  terminal  velocity  (lower  drag  force)  than  DELFIC 
for  the  larger  particles  (0.5  mm)  if  the  same  particle  density  is  used.  TASS,  however, 

assumes  a  lighter  density  for  the  dust  particles  (2.0  vs.  2.6  kg/m3)  which  compensates  for 
the  difference  in  the  drag  parameter.  The  stabilized  heights  for  the  same  particle  using  the 
two  codes  differ  by  less  than  100  meters,  therefore,  DELFIC’s  settling  subroutine  is  still 
used. 

Radial  Expansion.  The  last  major  modeling  effort  included  in  the  new  particle  rise 
subroutines  is  a  correction  for  the  radial  expansion  of  the  cloud.  All  of  the  velocity  equa¬ 
tions  presented  in  this  chapter  thus  far  are  for  a  vortex  sphere  of  constant  size  rising  in  the 
air.  Since  DELFIC’s  bubble  is  actually  expanding  continuously  during  its  rise,  a  correction 
is  needed  for  the  particle’s  velocity  equations. 
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The  expansion  correction  in  this  research  is  taken  to  be  proportional  to  the  location  of 
the  particle  and  the  current  bubble  expansion  rates.  From  the  bubble  rise  history,  the  base 
rate,  top  rate,  and  radial  rate  are  determined  for  the  time  of  interest.  These  values  are  used 
to  determine  the  vertical  and  radial  expansion  rates  for  the  particles  within  the  visible 
oblate  cloud.  A  particle’s  radial  velocity  will  increase  by  the  fraction  of  its  radial  distance 
to  the  radius  of  the  cloud,  multiplied  bv  the  radial  expiu.o.on  rate.  Likewise,  particles  in 
the  cloud  will  expand  vertically  commensurate  with  their  distance  from  the  cloud  center  to 
the  top  or  base. 

Integration  of  the  Particle  Equations.  Once  all  the  velocity  components  are  resolved 
from  the  underlying  forces  and  corrections,  the  particle’s  trajectory  is  calculated  by  inte¬ 
gration.  The  method  of  integration  chosen  is  a  fifth-order  Runge-Kutta  with  an  adaptive 
stepsize  to  ensure  accuracy  and  optimum  run  times  (52:701-715).  The  numerical  method 
chosen  was  well  documented  and  its  performance  checked  with  analytically  solvable 
equations.  A  main  calling  routine  for  the  vortex  particle  rise  module  was  created  (VPRM) 
as  well  as  a  subroutine  to  calculate  the  derivatives  of  the  particle’s  trajectory  for  integra¬ 
tion  (DERIVS). 

Comparison  Cases  and  Types  of  Comparisons 

The  vortex  particle  rise  model  as  described  above  was  then  used  to  simulate  a  number 
of  different  cases.  The  cases  used  include  the  Single  Burst  Nuclear  Cloud  Database 
(SBNCD)  and  the  1979  DELFIC  t; .  case  These  cases  were  used  io  investigate  the  effect 
of  using  vortex  flow  in  modeling  ,v  abilized  nuclea:  ,  *  ...  u.  it»  -  See  Table  7  on 
page  69).  A  number  of  diffeient  wa>s  of  displaying  the  •  -  t 
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Tbble  7:  Vortex  Particle  Rise  Comparison  Cases 


Case 

Yield  (kt) 

Scaled  Height  of  Burst  (ft/kt 

91A01 

300 

0 

91A02 

300 

50 

91A03 

300 

120 

91B01 

100 

0 

91B02 

100 

50 

91B03 

100 

120 

91F01 

800 

0 

91F02 

800 

50 

91H01 

20 

0 

1979  Test  Case 

50 

0 

Particle  Plots.  One  form  of  output  is  plots  of  trace  particles  as  a  function  of  time  for 
different  size  classes.  Although  these  were  not  used  for  comparison  with  hydrocode  data 
or  the  current  CRM,  they  do  help  show  the  effect  vortex  rise  has  on  the  trajectories  of  trace 
particles.  The  effects  of  radial  expansion,  particle  settling,  and  the  toroidal  flow  are  easily 
seen  by  these  plots. 

Density  Contours.  The  primary  means  of  comparing  the  vortex  particle  rise  results 
with  the  hydrocode  data  are  with  density  contour  plots.  These  plots  have  all  been  run  for 
the  SBNCD  cases  and  were  available  on  the  DNA  CRAY  at  Los  Alamos.  The  plots  arc 
presented  in  the  next  cl.  ‘er  and  give  both  a  radial  and  vertical  description  of  the  mass 
loading  of  the  cloud  at  various  times  from  one  minute  to  stabilization.  Plots  for  other  times 
and  other  variables  are  also  on  the  CRAY  (9). 

Grid  Size  and  Contour  Levels.  All  the  density  contours  shown  in  the  next  chapter 
were  generated  by  using  500  meter  gnd  cells  in  both  radial  and  vertical  dimensions  and 
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summing  the  mass  of  the  tracer  particles  which  lie  within  the  cell.  From  this  data,  contour 
lines  were  calculated  which  show  the  mass  density  in  10"X  g/cc  where  x  is  an  integer  from 
10  down  to  the  highest  density  present.  It  is  currently  thought  that  the  visible  cloud  has  a 

density  of  around  10'8  g/cc  (6).  The  plots  shown  in  the  next  chapter  (TASS  and  DELFIC) 
are  for  soil  loading  only.  TASS  plots  for  water  content  are  also  available  but  are  not  con¬ 
sidered  in  this  research. 

Mass  vs.  Altitude  Loading.  The  primary  purpose  in  comparing  vortex  particle  rise  to 
the  CRM’s  current  particle  rise  was  to  see  if  the  particles  stabilize  at  different  altitudes.  In 
particular,  would  the  larger  particles  attain  higher  altitudes  similar  to  preliminary  compar¬ 
isons  with  TASS.  For  this  reason,  mass  versus  altitude  was  chosen  as  the  means  for  com¬ 
paring  the  current  CRM  to  the  vortex  particle  rise  results  Another  reason  mass  versus 
altitude  was  chosen  was  that  the  CRM.  as  it’s  currently  setup,  has  no  radial  distribution. 

Fallout  Contours  and  Radial  Distribution.  The  final  comparison  made  was  a  prelimi¬ 
nary  review  of  the  implications  vortex  particle  rise  has  on  one  set  of  fallout  contours. 
Since  the  primary  purpose  for  upgrading  DELFIC  is  to  improve  the  DoD’s  fallout  predic¬ 
tion  capability,  the  fallout  contour  is  of  utmost  concern.  To  fully  see  how  a  radial  distribu¬ 
tion  affects  fallout  contours  requires  an  investigation  of  the  DTM  and  OPM  of  DELFIC. 
Since  this  was  not  within  the  scope  of  the  current  research,  conclusions  drawn  can  only  be 
preliminary  in  nature. 

Link  to  DELFIC’s  DTM.  To  produce  fallout  contours  using  DELFIC’s  DTM  and 
OPM  required  giving  up  the  radial  definition  of  the  stabilised  parucle  cloud.  The  distribu¬ 
tion  of  stabilized  particles  was  loaded  into  DELFIC’s  current  disk-like  structure  and  writ- 
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ten  to  tape.  The  sizes  of  the  disks  (both  vertical  and  radial)  were  chosen  to  be  similar  to 
disks  created  by  the  CRM,  with  the  only  difference  being  the  amount  of  mass  in  each  disk. 
The  DTM/OPM’s  method  of  assigning  mass  (activity)  in  the  radial  dimension  was  used 
for  the  generation  of  fallout  maps.  The  OPM  uses  a  bivariate  gaussian  distribution  which 
is  distorted  elliptically  (downwind  and  crosswind)  based  on  diffusion  during  the  transport. 
The  difference  in  vertical  loading,  therefore,  is  the  only  difference  investigated  in  the  fall¬ 
out  contour  comparison  study. 


This  chapter  presents  the  results  and  comparisons  of  the  two  main  areas  of  research  \ 

conducted  by  the  author.  The  chapter  is  organized  into  two  sections,  bubble  rise  and  parti¬ 
cle  rise.  The  first  section  gives  the  results  of  comparing  the  bubble  rise  subroutines  to 
observation.  The  second  main  section  gives  a  comparison  of  vortex  theory  particle  rise  to 
the  current  CRM  particle  rise  and  to  hydrocode  results. 

Bubble  Rise 

This  section  presents  the  results  obtained  from  the  different  versions  of  the  CRM 
described  in  “Work  Done  on  DELFIC’s  Bubble  Rise  Subroutines”  on  page  31 .  The  fol¬ 
lowing  subsections  show  the  performance  of  the  1979,  corrected,  and  improved  versions 
of  the  CRM.  The  1979  version  is  the  documented  code  as  provided  to  the  author.  The  cor¬ 
rected  version  is  the  code  with  the  programming  errors  removed.  The  improved  version  is 
the  code  with  the  single  term  entrainment  equation  and  constant  parameters.  The  three  dif¬ 
ferent  versions  are  compared  to  actual  observations  from  the  US  atmospheric  nuclear  test 
events.  This  is  followed  by  a  discussion  of  observation  error  and  a  discussion  of  the  figure 
of  merit  chosen.  Finally,  a  couple  of  examples  of  the  time-dependent  performance  of  the 
improved  CRM  are  given. 

The  performance  of  stabilized  bubble  prediction  lies  in  matching  observed  nuclear 
cloud  values,  in  particular  stabilized  cloud  tops.  Norment  points  out  that  the  cloud  top.  by 
its  very  nature,  is  crucial  to  the  art  of  predicting  fallout,  the  ultimate  goal  of  DELFIC  The 
cloud  top  value  determines  what  wind  layer  values  are  to  be  used  in  transporting  the  radio¬ 
activity.  It  is  also  noted  that  cloud  top  recordings  are  by  far  the  most  reliable  data  from  the 
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atmospheric  nuclear  tests.  Cloud  bases  and  radii,  on  the  other  hand,  are  ill-defined  due  to 
the  presence  t  f  a  stem  and  the  wind  shearing  of  the  clouds.  Norment  emphasizes  the  fact 
that  the  least  important  property  is  the  observed  cloud  base,  since  stabilized  particle  posi¬ 
tions  in  the  cloud  may  not  correspond  with  this  visible  base  value  (46:21).  Of  course  tem¬ 
poral  performance  of  the  rising  cloud  needs  to  be  considered  and  is  also  reported  here. 

The  decision  was  made  to  repeat  the  validation  of  the  1979  and  newer  versions  of  the 
CRM  in  the  same  manner  as  Nonnent’s  1977  validation  effort.  Therefore,  the  figures  of 
merit  used  were  the  fractional  root  mean  square  erTor  (FRMS),  and  the  fractional  mean 
deviation  (FMD).  Both  were  described  in  the  section  titled  “Fractional  Root  Mean  Square 
(FRMS)  Error  Definition”  on  page  38.  It  should  be  noted  at  this  point  that  the  FMD  for  the 
entire  set  of  validation  shots  may  be  deceiving  since  overpredictions  can  balance  out 
undeipredictions.  The  deviation  for  an  individual  shot  is  of  value  since  it  tells  if  the  model 
overpredicts  or  underpredicts  the  height,  depending  on  its  sign  (negative  and  positive 
respectively). 

The  FRMS  is  chosen  as  the  primary  figure  of  merit  since  it  sums  absolute  error  in 
either  direction.  The  combination  of  the  FRMS  applied  to  the  cloud  top  values  was  chosen 
to  adjust  the  constants  in  the  improved  model  (to  be  discussed  later).  The  results  of  pre¬ 
dicting  cloud  bases  are  shown  since  observation  values  were  available  for  most  of  the 
shots  (49  of  54).  Radii  observations  were  not  widely  available,  and  so  are  not  presented  at 
all.  Table  8  on  page  74  and  Table  9  on  page  76  present  the  results  of  the  bubble  rise  part  of 
the  research.  The  values  shown  are  for  the  1979  version  of  the  CRM  as  well  as  the  cor¬ 
rected  and  improved  versions  of  the  code.  They  are  all  relative  to  the  burst  height  (HOB 
plus  surface  altitude). 
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Table  8:  Cloud  Tbp  Comparison  of  Models  to  Observation  (relative  to  burst  point) 


Shot 

Yield 

(kt) 

Observed 
Cloud 
Tbp  (m) 

Cal  cl 

1979 

lated  Clou 

Corrected 

d  Tbp  (m) 

Improved 

Fr 

1979 

actional  Di 

Corrected 

sviation 

Improved 

Humboldt 

0.0078 

1050 

759 

827 

1038 

0.28 

0.01 

Catron 

SHI 

1344 

1061 

1170 

1415 

0.21 

-0.05 

Vesta 

1760 

2282 

2294 

2247 

-0.30 

-028 

DonaAna 

0.037 

1940 

2831 

2802 

2796 

-0.46 

-0.44 

-0.44 

Hidalgo 

0.077 

2267 

2172 

2258 

2515 

0.04 

0.00 

-0.11 

Quay 

0.079 

1722 

1548 

1543 

1768 

0.10 

0.10 

-0.03 

Eddy 

0.083 

1925 

2368 

2137 

2635 

-0.23 

-0.11 

-0.37 

RioArriba 

0.09 

2870 

1895 

1975 

2559 

0.34 

0.31 

0.11 

Wrangell 

0.115 

1653 

1627 

1640 

1861 

0.02 

0.01 

-0.13 

Franklin 

0.14 

3772 

4175 

4253 

4179 

-0.11 

-0.13 

-0.11 

Wheeler 

0.197 

3740 

3374 

3396 

3684 

0.10 

0.09 

0.01 

Ray 

2644 

2121 

2141 

2524 

0.20 

0.19 

0.05 

Ruth 

0.2 

2833 

2838 

2946 

3051 

0.00 

-0.04 

-0.08 

JohnnieBoy 

0.5 

3612 

2575 

2565 

2818 

0.29 

0.29 

0.22 

Laplace 

1 

4592 

4709 

4736 

4819 

-0.03 

-0.03 

-0.05 

SantaFe 

1.3 

3753 

3550 

3569 

4254 

0.05 

0.05 

-0.13 

Lea 

1.4 

3449 

3931 

3925 

4536 

-0.14 

-0.14 

-0.32 

Mora 

3906 

4339 

4374 

4521 

-0.11 

-0.12 

-0.16 

John 

6008 

4197 

4197 

4540 

0.30 

0.30 

0.24 

DeBaca 

3601 

3878 

3811 

4989 

-0.08 

-0.06 

-0.39 

FranklinPrime 

8249 

5467 

5480 

5841 

0.34 

0.34 

0.29 

Sanford 

4.9 

6530 

4946 

4986 

5942 

0.24 

0.24 

0.09 

Socorro 

6 

6207 

5776 

5792 

6287 

0.07 

0.07 

-0.01 

Morgan 

8 

10755 

6374 

6379 

6773 

0.41 

0.41 

0.37 

Owens 

9.7 

9231 

8033 

7999 

7407 

0.13 

0.13 

0.20 

Wilson 

10 

9226 

6409 

6389 

7987 

0.31 

0.31 

0.13 

Kepler 

10 

7069 

7612 

7611 

7648 

-0.08 

-0.08 

-0.08 

Fizeac 

11 

10811 

7780 

7833 

7915 

0.28 

0.28 

0.27 

liable  8:  Cloud  lbp  Comparison  of  Models  to  Observation  (relative  to  burst  point) 


Shot 

Yield 

(kt) 

Observed 
Cloud 
Tbp  (m) 

Calcu 

1979 

lated  Clou 

Corrected 

d  Tbp  (m) 

Improved 

Fr 

1979 

actional  Dc 

Corrected 

aviation 

Improved 

Galileo 

11 

9830 

7554 

7573 

8029 

0.23 

0.18 

Doppler 

11 

9836 

7731 

7685 

7307 

0.22 

0.26 

Dixie 

11 

10654 

8099 

8092 

8681 

0.24 

0.24 

0.19 

Boltzman 

12 

8615 

10524 

10517 

9887 

-0.22 

-0.22 

-0.15 

Newton 

12 

8021 

7958 

8008 

8057 

0.01 

000 

0.00 

Charleston 

12 

8012 

6779 

6823 

6801 

0.15 

0.15 

0.15 

Grable 

15 

9570 

6164 

6147 

6425 

0.36 

0.36 

0.33 

Annie 

16 

11178 

9365 

9480 

10039 

0.16 

0.15 

0.10 

Shasta 

17 

8264 

9347 

9215 

8717 

-0.13 

-0.12 

-0.05 

Diablo 

17 

8239 

8884 

8837 

9171 

-0.08 

-0.07 

-0.11 

Whitney 

19 

7624 

8459 

8474 

9042 

-011 

-0.11 

-0.19 

Stokes 

19 

9545 

8494 

8497 

8732 

0.11 

0.11 

0.09 

Badger 

23 

9513 

7554 

7568 

8897 

0.21 

0.20 

0.06 

Nancy 

24 

11244 

8807 

8823 

9217 

0.22 

0.22 

0.18 

Encore 

27 

11125 

8908 

8911 

9246 

0.20 

0.20 

0.17 

Hanry 

32 

11642 

11844 

11997 

12640 

-0.02 

-003 

-0.09 

Priscilla 

37 

11955 

10782 

10774 

11150 

0.10 

0.10 

0.07 

Lacrosse 

40 

11582 

7410 

9164 

8983 

0.36 

0.21 

0.22 

Simon 

43 

12028 

12140 

Smoky 

44 

10004 

11290 

Climax 

61 

11382 

12084 

Hood 

74 

12884 

12719 

Koon 

110 

16150 

15549 

Zuni 

3500 

24076 

25195 

Tewa 

5000 

30171 

26613 

Bravo 

15000 

34745 

35450 

12183  12181  -0.01  -0.01  -0.01 

11298  11181  -0.13  -0.13  -0.12 

12092  12053-0.06  -0.06  -0.06 

12724  13245  0.01  0.01  -0.03 

15713  14995  0.04  0.03  0.07 

27341  27282  -0.05  -0.14  -0.13 

29399  29520  0.12  0.03  0.02 

37085  36118  -0.02  -0.07  -0.04 

FMD  0.08  0.06  0.01 

FRMS  0.20  0.19  0.18 


Table  9:  Cloud  Base  Comparison  of  Models  to  Observation  (relative  to  burst  point) 


Shot 

Yield 

<kt) 

Observed 

Cloud 

Base(m) 

Caicul 

1979 

lated  Cloud 

Corrected 

1  Base  (m) 

Improved 

Fr 

1979 

actional  De 

Corrected 

aviation 

Improved 

Humboldt 

593 

450 

471 

688 

Catron 

0.021 

277 

646 

724 

957 

.  1 

-2.45 

Vesta 

1577 

1601 

1642 

DonaAna 

0.037 

568 

1657 

1682 

1942 

•1.92 

-2.42 

Hidalgo 

0.077 

1048 

1297 

1381 

1740 

-0.24 

-0.66 

Quay 

0.079 

960 

921 

923 

1219 

0.04 

0.04 

-0.27 

Eddy 

0.083 

858 

1477 

1419 

1634 

-0.72 

-0.65 

-0.90 

RioArriba 

0.09 

2108 

1187 

1238 

1699 

0.44 

0.41 

0.19 

Wrangell 

0.115 

739 

984 

994 

1289 

-0.33 

-0.35 

-0.75 

Franklin 

0.14 

2949 

2674 

2804 

2984 

0.09 

0.05 

-0.01 

Wheeler 

0.197 

2825 

2104 

2114 

2539 

0.26 

0.25 

0.10 

Ray 

0.2 

1089 

1326 

1343 

1676 

-0.22 

-0.23 

-0.54 

Ruth 

0.2 

1949 

1551 

1705 

2147 

0.20 

0.13 

-0.10 

JohnnieBoy 

0.5 

2240 

1629 

1638 

1921 

0.27 

0.27 

0.14 

Laplace 

1 

2763 

3054 

3077 

3417 

-0.11 

-0.11 

-0.24 

SantaFe 

1.3 

2229 

2230 

2226 

2899 

0.00 

0.00 

-0.30 

Lea 

1.4 

1925 

2468 

2468 

3054 

-0.28 

-0.28 

-0.59 

Mora 

2 

1315 

2845 

2871 

3150 

-1.16 

-1.18 

-1.40 

John 

2 

2702 

2708 

3100 

DeBaca 

2.2 

1315 

2372 

2358 

3367 

-0.80 

-0.79 

-1.56 

FranklinPrime 

4.7 

4896 

3601 

3605 

3979 

0.26 

0.26 

0.19 

Sanford 

4.9 

2415 

3112 

3129 

3827 

-0.29 

-0.30 

-0.58 

Socorro 

6 

4378 

3825 

3839 

4292 

0.13 

0.12 

0.02 

Morgan 

8 

6488 

4183 

4187 

4594 

0.36 

0.35 

0.29 

Owens 

4659 

5464 

5442 

5107 

-0.17 

-0.17 

-0.10 

Wilson 

■ 

6178 

4209 

4193 

5169 

0.32 

0.32 

0.16 

Kepler 

■ 

46:0 

5100 

5077 

5172 

-0.10 

-0.10 

-0.12 

Fizeau 

■ 

6849 

5241 

5268 

5451 

0.23 

0.23 

0.20 
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liable  9:  Cloud  Base  Comparison  of  Models  to  Observation  (relative  to  burst  point) 


Shot 

Yield 

<kt) 

Observed 

Cloud 

Base(m) 

Calcu 

1979 

lated  Clouc 

Corrected 

1  Base  (m) 

Improved 

Ft 

1979 

Galileo 

11 

3734 

5021 

5471 

-0.34 

-0.47 

Doppler 

11 

5264 

5027 

4940 

0.04 

2 

0.06 

Dixie 

11 

6996 

5251 

5713 

0.25 

0.25 

0.18 

Boltzman 

12 

5567 

7258 

7246 

6759 

-0.30 

-0.30 

-0.21 

Newton 

12 

4058 

5181 

5218 

5411 

-0.28 

-0.29 

-0.33 

Charleston 

12 

4354 

4643 

4684 

4806 

-0.07 

-0.08 

-0.10 

Grable 

15 

5913 

3866 

3861 

4233 

0.35 

0.35 

0.28 

Annie 

16 

7216 

5917 

5953 

6614 

0.18 

0.18 

0.08 

Shasta 

17 

3387 

6118 

6036 

5872 

-0.81 

-0.78 

-0.73 

Diablo 

17 

4581 

5960 

5920 

6306 

-0.30 

-0.29 

-0.38 

Whitney 

19 

3967 

5570 

5568 

6133 

-0.40 

-0.40 

-0.55 

Stokes 

19 

6497 

5449 

5439 

5758 

0.16 

1 

0.16 

0.11 

Badger 

23 

5550 

4758 

4772 

5762 

0.14 

0.14 

-0.04 

Nancy 

24 

6520 

5782 

5801 

6214 

0.11 

0.11 

0.05 

Encore 

27 

6858 

5474 

5476 

5981 

0.20 

0.20 

0.13 

Harry 

32 

7070 

7213 

7377 

8105 

-0.02 

-0.04 

-0.15 

Piiscilla 

37 

6164 

6914 

6909 

7392 

-0.12 

-0.12 

-0.20 

Lacrosse 

40 

6096 

5180 

5856 

6011 

0.15 

0.04 

0.01 

Simon 

43 

8065 

7907 

7948 

8090 

0.02 

0.01 

0.00 

Smoky 

44 

7566 

7566 

7677 

Climax 

61 

9035 

7883 

7891 

8141 

0.13 

0.13 

0.10 

Hood 

74 

8921 

8065 

8065 

8737 

0.10 

0.10 

0.02 

Koon 

110 

9875 

10035 

9649 

Zuni 

3500 

14932 

13623 

15016 

15374 

0.09 

-0.01 

-0.03 

Tewa 

5000 

14017 

16369 

) 

16435 

Bravo 

15000 

16853 

20938 

22055 

20240 

-0.24 

-0.31 

-0.20 

FMD 

-0.12 

-0.14 

-0.29 

FRMS 

0.47 

0.49 

0.66 
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Observation  The  observed  values  of  the  cloud  tops  were  taken  from  Refs  16  and  17, 
and  for  the  majority  were  given  to  the  nearest  kilofoot  (305  m).  The  values  for  cloud  tops 
and  bases  are  plotted  as  a  function  of  yield  in  Figure  13  on  page  79.  No  mention  was  made 
to  the  estimated  error  in  the  observation  in  these  references.  A  search  was  made  of  the  old 
weapon  test  (WT)  reports  to  try  and  get  a  rough  idea  of  the  appropriate  error,  but  to  no 
avail. 

One  study  did  comment  on  reasons  why  early  data  was  considered  questionable.  It 
stated  that  some  data  were  considered  unreliable  due  to  the  lack  of  a  drift  correction  when 
the  cloud  vectored  toward  the  camera.  Some  data  were  inaccurate  due  to  manual  operation 
of  the  theodolite  equipment  or  conflicting  results  with  the  camera  data  (14: 14). 

To  get  a  feeling  for  the  impact  of  observation  error  on  the  FRMS  value,  a  simple  sensi¬ 
tivity  study  was  done.  For  the  purpose  of  this  side  study,  it  was  assumed  that  all  shots  in 
the  validation  set  contained  the  same  amount  of  absolute  error  and  the  FRMS  was  calcu¬ 
lated.  The  FRMS  associated  with  this  exact  amount  of  error  is  shown  along  with  the  value 
of  the  FMD  for  the  entire  set  in  Table  10  on  page  80.  As  can  be  seen  from  the  table,  if  one 
assumes  the  observed  values  can  be  off  by  one  kilofoot  (305  m),  the  FRMS  is  already 
approaching  0.1. 

Norment  presented  results  of  his  validation  study  using  the  FRMS  and  FMD  values 
(46:27).  These  numbers  are  shown  in  Table  11  on  page  80.  Included  in  the  table  are  the 
values  for  a  set  of  53  shots  (54  minus  Humboldt)  and  60  shots  (54  plus  6  classified  shots). 
Also  listed  are  the  results  obtained  when  the  author  ran  the  tabulated  values  in  Norment’s 
1970  simulation  study  of  the  54  shots.  Norment  did  not  calculate  his  own  value  for  FRMS 
with  the  1970  version  of  the  CRM. 
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Table  10:  Sensitivity  of  FRMS  and  FMD  to  a  Constant  Error  in  the  Cloud  Top  Value 


Assumed  Error  (ft) 

FRMS 

FMD 

1000 

0.089 

0.066 

2000 

0.177 

0.132 

3000 

0.266 

0.197 

4000 

0.354 

0.263 

5000 

0.443 

0.329 

Table  11:  Cloud  Top  FRMS  and  FMD  for  the  1970  &  1977  Versions  of  the  CRM 


1970  (54  shots) 

1977  (53  shots) 

1977  (60  shots) 

FMD 

0.16 

-0.009 

0.002 

FRMS 

0.24 

0.15 

0.14 

1979  Version.  First  the  1979  version  of  the  CRM  was  used  to  calculate  the  cloud  tops 
and  bases  of  the  validation  set.  The  results  of  these  runs  are  plotted  in  Figure  14  on 
page  81  as  well  as  being  listed  in  Tabic  8  on  page  74  and  Table  9  on  page  76  (which  show 
the  values  for  the  FRMS  figure  of  merit). 

The  first  question  that  comes  to  mind  is  why  the  1979  version  of  the  model  performs 
worse  than  the  previously  tabulated  values  for  the  1977  version.  Since  Norment  did  not 
tabulate  the  data  used  to  compute  the  reported  FRMS  values  in  1977,  there  is  no  way  of 
knowing  why  the  1979  CRM  performs  worse.  The  1979  CRM  was  checked  line  by  line  by 
the  author  and  the  test  case  cloud  top  listed  in  the  documentation  was  duplicated  tc  all  four 
significant  digits.  Therefore,  the  1979  FRMS  value  is  used  as  a  baseline  for  compaiison, 
and  the  1977  values  are  held  suspect. 

Corrected  Version.  After  making  the  corrections  mentioned  in  “Corrections  Needed 
in  the  1979  DELFIC  CRM”  on  page  31,  the  set  of  54  shots  was  again  run.  The  results  of 
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Figure  14:  1979  CRM  Cloud  Tops  and  Bases  (relative  to  burst  heights) 
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this  newer  version  of  the  CRM  are  plotted  in  Figure  15  on  page  83.  As  can  be  seen  from 
Table  8  on  page  74,  the  FRMS  value  decreases  slightly  from  0.20  to  0. 19  after  making  the 
suggested  corrections. 

Improved  Version.  The  improved  version  of  the  CRM  represents  a  change  back  to  a 
single  term  entrainment  equation  and  constant  values  for  the  entrainment  (|i)  and  eddy  vis¬ 
cous  drag  (k2)  parameters.  The  optimum  values  for  the  parameters  were  found  by  running 
a  matrix  of  values  for  the  pair.  The  suggested  ranges  by  Huebsch  (21:72)  were  tried  ini¬ 
tially.  His  range  for  the  entrainment  parameter  (0. 15  to  0.25)  was  too  high,  most  likely  due 
to  the  model’s  change  in  cloud  shape  modeling  (See  “Cloud  Shape  and  Volume”  on 
page  30).  A  range  of  smaller  values  (0.08  to  0.15)  closer  to  the  values  used  in  the  1979 
version  was  tried.  As  for  the  eddy  viscous  drag  parameter,  values  from  Huebsch ’s  highest 
recommended  value  (0.125)  down  to  0.0  were  tried. 

The  FRMS  results  of  varying  the  entrainment  and  eddy  viscous  drag  parameters  are 
shown  in  Figure  16  on  page  84.  As  can  be  seen  from  the  saddle  shape  of  the  figure,  there  is 
a  compensating  effect  when  lowering  the  value  of  one  parameter  while  raising  the  other.  It 
was  found  that  the  minimum  value  for  the  ranges  considered  was  FRMS=0.18.  This  value 
is  achieved  with  more  than  one  combination  of  parameters.  The  pair  chosen  was  p=0.12 
and  £2=0. 1 .  These  values  correspond  to  the  1979  CRM  values  for  yields  up  to  a  kiloton 
(See  Figure  4  on  page  28). 


The  results  of  using  these  values  along  with  a  single  term  entrainment  equation  are 
shown  in  Figure  17  on  page  85.  As  can  be  seen  from  Table  8  on  page  74,  the  FRMS  value 
decreases  again  from  0.19  tc  0.18  after  making  these  improvements.  This  stands  in  con- 
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Figure  16:  FRMS  Values  for  Entrainment  (p.)  and  Eddy  Viscous  Drag  (£3)  Parameter 

Combinations 


trast  to  comments  made  in  Norment’s  1977  report:  “Extensive  calculations... have  shown 
that  it  is  not  possible  to  adequately  »..atch  calculated  with  observed  results  if  single,  yield 
independent  values  are  used  for  these  [(!  and  parameters”  (46:40).  One  could  argue 
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that  the  FRMS  with  constant  parameters  (See  Figure  16  on  page  84)  is  still  not  "ade¬ 
quate”,  but  it  is  certainly  an  improvement. 

The  FMD  also  reduces  substantially  for  the  whole  group.  However  this  may  imply  that 
there  are  just  as  many  undeipredictions  as  overpredictions.  But  this  near  zero  FMD  com¬ 
bined  with  the  fact  that  the  FRMS  is  reduced,  implies  a  somewhat  better  fit  to  observed 
cloud  tops  with  the  chosen  yield  independent  parameters.  The  figure  of  merit  for  the  cloud 
bases  on  the  other  hand  gets  about  30-40%  worse.  Considering  the  scatter  in  the  observed 
data  (See  Figure  13  on  page  79),  and  the  fact  that  cloud  bases  are  ill-defined,  not  too  much 
can  be  implied  from  this  result. 

Time  Dependent  Results.  The  temporal  performance  of  the  bubble  rise  is  the  final 
area  to  be  discussed.  Although  the  figures  of  merit  chosen  were  not  related  to  the  rate  of 
rise,  the  way  a  cloud  gets  to  its  final  height  does  play  a  part  in  determining  the  strength  of 
the  flow  in  the  particle  rise  modeling.  Figure  18  on  page  87  and  Figure  19  on  page  88 
show  the  calculated  time  history  plots  of  shots  Castle  Bravo  and  Upshot-Knothole  Simon. 
These  cases  were  chosen  since  Norment  used  them  in  his  1977  study  and  temporal  plots 
were  available  from  DASA  1251-EX  (16;  17).  These  references  contain  temporal  plots  for 
16  of  the  54  shots  in  the  set. 

While  the  results  for  both  Simon  and  Bravo  are  excellent  as  far  as  matching  cloud 
tops,  both  achieve  the  desired  value  earlier  than  observation.  The  Bravo  simulation  (See 
Figure  18  on  page  87)  shows  the  cloud  top  stabilizing  at  3  minutes,  while  the  observation 
(See  Figure  20  on  page  89)  puts  vertical  stabilization  at  4.5  minutes.  The  observation  also 
shows  oscillation  in  the  cloud  top  height  following  the  apogee,  which  is  typical  for  high 
yield  thots.  The  observed  cloud  base  continues  to  rise  even  at  10  minutes  while,  by  the 
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Figure  18:  Calculated  Cloud  Dimensions:  Operation  Castle  Shot  Bravo 


nature  of  the  CRM  modeling,  the  simulated  cloud  base  is  made  to  stop  rising  at  the  same 
time  as  the  top. 

The  simulation  of  Simon  (See  Figure  19  on  page  88)  shows  the  cloud  top  achieving 
vertical  stabilization  at  7  minutes  while  observation  (See  Figure  21  on  page  90)  puts  verti¬ 
cal  stabilization  at  10  minutes.  With  Simon,  however,  the  observed  base  vertically  stabi¬ 
lizes  at  7  minutes,  identical  to  the  simulation.  It  should  be  noted  that  the  simulated  vertical 
stabilization  times  for  both  of  these  shots  was  about  the  same  as  those  simulated  in  Nor- 
ment’s  1977  study.  While  these  two  simulations  both  stabilize  the  cloud  top  earlier  than 
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Figure  19:  Calculated  Cloud  Dimensions:  Operation  Upshot-Knothole  Shot  Simon 


observation,  other  simulations  (e  g.  Upshot  Knothole  events  Annie  and  Ruth)  stabilize  the 
cloud  top  later  than  observation. 

Particle-Rise 

This  second  major  section  of  the  chapter  presents  the  results  of  running  the  particle 
rise  subroutines  which  use  the  vortex  flow  field  equations.  There  are  two  main  compari¬ 
sons  to  be  shown:  1)  vortex  particle  rise  to  hydrocode  dust  modeling,  and  2)  vortex  parti¬ 
cle  rise  to  DELFlC's  1-D  particle  rise.  Before  these  comparisons  are  made,  however, 
sample  snapshots  of  the  trace  particles  ’  histories  are  given  for  physical  insight. 
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Figure  20:  Observed  Cloud  Dimensions:  Operation  Castle  Shot  Bravo  (17:68) 


Particle  Plots.  Animation  of  the  trace  particles  is  a  useful  tool  for  evaluating  the 
code’s  functioning  and  overall  performance.  Animation  in  Mathematica  was  used  as  an 
integral  part  of  the  particle  rise  research.  It  allowed  the  author  to  view  the  effects  of  the 
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Figure  21:  Observed  Cloud  Dimensions:  Operation  Upshot-Knothole  Simon  (16:141) 


submodel  choices  such  as  gravity  settling,  radial  expansion,  and  initial  source  location. 


Although  animation  is  not  able  to  be  placed  in  this  report,  a  sample  of  a  sequence  of  snap¬ 
shots  can  be  given. 

In  Figure  22  on  page  92  the  positions  of  trace  particles  are  shown  from  initial  time 
(4.35  sec)  to  stabilization  time  (608.1  sec)  along  with  the  whole-minute  snapshots.  The 
sequence  represents  the  1979  test  case  using  the  improved  CRM.  For  this  figure  there  are 
5  particle  size  groups  (22, 60, 130,  277,  and  768  microns),  each  with  a  20  by  20  grid  of 
starting  locations,  for  a  total  of  2000  trace  particles.  One  must  keep  in  mind  that  the  trace 
particles  do  not  all  have  equal  weights  (of  mass  or  radioactivity).  This  is  due  to  the  fact 
that  each  trace  particle  represents  an  annulus  at  the  radius  where  it  originates.  Therefore, 
particles  starting  at  the  outer  radii  have  a  larger  weight  than  those  starting  at  the  inner 
radii. 

The  first  few  snapshots  show  how  the  vortex  flow  quickly  moves  the  dust  out  into  a 
torus  away  from  the  axis.  In  addition,  one  can  see  the  number  of  trace  particles  in  the  cen¬ 
ter  of  the  torus  is  sparser  than  at  the  surface  of  the  torus.  In  the  second  row  of  snapshots, 
the  reader  can  see  that  the  cloud  has  reached  vertical  stabilization.  After  this,  radial  expan¬ 
sion  is  seen,  as  the  particle  sizes  also  start  to  separate.  The  larger  particles  fall  first  and 
some  exit  the  cloud  before  stabilization.  Once  they  fall  below  the  cloud  base  they  no 
longer  are  affected  by  the  expansion  correction.  The  last  snapshot  shows  the  different 
groups,  with  the  two  smallest  groups  overlapping  each  other.  This  occurs  because  neither 
of  these  smaller  groups  has  an  appreciable  terminal  velocity. 

Dust  Density  Contour  Plots.  The  first  set  of  comparisons  to  be  made  is  to  the  hydro- 
code  dust  plots.  Each  of  the  9  different  cases  described  in  Table  7  on  page  69  are  plotted. 
Plots  are  provided  for  both  DELFIC  with  vortex  flow  and  TASS  hydrocode  results.  Each 
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Figure  22:  Cross  Section  of  Particle  Cloud  with  5  Trace  Particle  Sizes  as  it  Rises 


(snapshots  every  minute:  initial  time  to  stabilization  time;  distances  in  kilometers) 
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required  multiple  plots  per  simulation  for  times  from  one  to  ten  minutes.  The  plots  are 
shown  in  Figure  23  on  page  94  through  Figure  40  on  page  1 11  for  comparison.  Shown  are 
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the  total  dust  contours  at  levels  of  10*l°  g/cc  and  higher  by  factors  of  10.  These  simula¬ 
tions  were  run  with  the  full  number  of  particle  size  groups  shown  in  Table  6  on  page  43. 

DELFIC  with  Vortex.  The  radioactive  particles  being  simulated  quickly  rise  up 
and  out  into  a  torus.  There  is  only  a  small  amount  of  this  dust  which  shows  up  in  the  stem 
or  pedestal  at  early  times.  After  vertical  stabilization,  radial  expansion  becomes  the  domi¬ 
nant  force  along  with  gravity  settling.  The  dust  contours  above  the  main  cloud  base 
expand  out,  emphasizing  the  void  on  the  axis  of  the  cloud.  Particle  settling  of  the  larger 
groups  does  form  a  visible  stem  at  late  times. 

These  plots  show  that  the  dust  that  originates  in  the  cloud  at  4-5  sec,  is  not  uniformly 
distributed  radially  at  later  times.  When  the  cloud  stops  rising,  the  vortex  flow  stops  and 
the  particles  move  under  the  influence  of  gravity  and  expansion  only.  This  phase  of  the 
particle  modeling  is  limited  in  that  it  does  not  account  for  any  diffusive  motion  of  the  par¬ 
ticles  back  to  the  areas  of  lower  concentration.  It  also  turns  off  the  toroidal  flow  without 
allowing  for  a  residual  continuation  of  this  flow  after  vertical  stabilization. 

TASS  Output.  The  hydrocode  output  does  account  for  both  diffusion  and  a  grad¬ 
ual  decay  of  the  flow  field  velocity  after  vertical  stabilization.  Tliis  residual  motion  may 
account  for  a  retardation  of  vertical  settling  in  the  TASS  results.  Also  noticeably  different 
is  the  presence  of  an  appreciable  stem  throughout  the  history  of  the  dust  cloud.  The  TASS 
dust  clouds  also  reach  a  higher  altitude  than  DELFIC,  especially  the  highest  near-suiface 
bursts  with  SHCB=120  ft/kt1^  (See  “Appendix  A:  Cloud  Tops”  on  page  127). 
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Figure  23:  91H01  (20  kt  0  SHOB)  DELFIC  Density  Contours  (10*  g/cc) 
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Figure  24:  91H01  (20  kt  0  SHOB)  TASS  Density  Contours  (10*  glee) 
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Figure  25:  91B01  (100  kt  0  SHOB)  DELF1C  Density  Contours  (10'*  g/cc) 
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Figure  26:  91B01  (100  kt  0  SHOB)  TASS  Density  Contours  (10*  g/cc) 
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Figure  27:  91B02  (100  kt  50  SHOB)  DELFIC  Density  Contours  (10*  g/cc) 
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Figure  28:  91B02  (100  kt  50  SHOB)  TASS  Density  Contours  (10‘x  g/cc) 
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Figure  29:  91B03  (100  kt  120  SHOB)  DELFIC  Density  Contours  (10*x  g/cc) 
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Figure  30:  91B03  (100  kt  120  SHOB)  TASS  Density  Contours  (10'x  g/cc) 
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Figure  31:  91A01  (300  kt  0  SHOB)  DELFIC  Density  Contours  (10'x  g/cc) 
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Figure  32:  91A01  (300  kt  0  SHOB)  TASS  Density  Contours  (10'x  g/cc) 
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Figure  33:  91A02  (300  kt  50  SHOB)  DELFIC  Density  Contours  (10'*  glee ) 
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Figure  34:  91A02  (300  kt  50  SHOB)  TASS  Density  Contours  (10'x  g/cc) 
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Figure  35:  91A03  (300  kt  120  SHOB)  DELFIC  Density  Contours  (10*x  g/cc) 
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Figure  36:  91A03  (300  kt  120  SHOB)  TASS  Density  Contours  (10  *  g/cc) 
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Figure  37:  91F01  (800  kt  0  SHOB)  DELFIC  Density  Contours  (10'x  g/cc) 
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Figure  38:  91F01  (800  kt  0  SHOB)  TASS  Density  Contours  (10*  g/cc) 


109 


DST+PBL 


91F02  DELFIC 


0.0  2.5 


DST+PBL 


2.5  s.o  9:8  2-5  5*° 

RADIUS  (KM) 

91F02  DELFIC 


2.5  5.0  7.5  iO.O 


RADIUS  (KM) 


Figure  39:  91F02  (800  kt  50  SHOB)  DELFIC  Density  Contours  (10x  g/cc) 
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Figure  40:  91F02  (800  kt  50  SHOB)  TASS  Density  Contours  (10'x  g/cc) 
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A  key  difference  again  is  that  the  TASS  dust  clouds  represent  both  radioactive  dust 
and  nonradioactive  dust  alike.  The  hydrocodes  also  account  for  a  time  dependent,  empiri¬ 
cal  loading  of  the  dust  cloud  from  two  sources,  ejecta  and  sweepup.  The  source  of  dust 
from  the  surface  which  originates  at  late  times  after  the  cloud  has  lifted  off  the  ground, 
(i.e.  sweepup),  may  account  for  the  dust  in  the  stem  and  close  to  the  axis  of  the  cloud.  This 
dust  may  not  enter  the  torus  due  to  a  reduced  vortical  intensity  at  late  times.  A  partitioning 
of  the  radioactive  dust  within  the  main  cloud  may  then  result  with  the  torus  containing 
almost  all  of  the  radioactivity.  There  currently  is  no  way  of  distinguishing  radioactive  dust 
from  nonradioactive  dust  in  the  MAZ/TASS  hydrocode  results. 

Vertical  Distribution  of  Mass.  The  second  set  of  comparisons  to  be  made  in  this  sec¬ 
tion  is  a  comparison  of  the  1-D  particle  rise  to  the  vortex  flow  method.  The  primary  focus 
will  be  on  the  vertical  distribution  of  different  particle  sizes.  The  radial  distribution  will 
also  be  discussed  with  implications  on  fallout  contours  mentioned. 

Plots  of  mass  as  a  function  of  altitude  are  shown  for  the  1979  test  case  (See  Figure  41 
on  page  1 13).  Results  from  both  the  old  and  new  methods  of  particle  rise  are  shown  in 
each  frame  of  the  figure  (DELFIC’s  1-D  on  the  left  and  vortex  theory  on  the  right).  Values 
shown  are  the  fractions  of  the  total  mass  for  the  vertical  height  of  the  cloud  layer  shown  in 
kilometers.  The  horizontal  lines  represent  the  CRM  simulated  cloud  top  and  base  alti¬ 
tudes.  The  first  frame  represents  all  sizes  in  the  total  cloud  while  the  other  frames  are  for 
specific  particle  sizes. 

One-Dimensional  Results.  The  CRM  currently  assumes  the  cloud  starts  off  uni¬ 
formly  loaded  at  the  beginning  of  its  rise.  The  smaller  groups,  which  have  no  appreciable 
settling  velocity,  are  therefore  still  close  to  uniform  in  their  vertical  distribution  at  stabili- 


112 


Figure  41:  Stabilized  Cloud  Mass- Altitude  Loading  for  1979  Test  Case  (without 
vortex  modeling  on  left  and  with  vortex  modeling  on  right  for  5  particle  sizes) 


zation.  The  heavier  groups  do  have  a  significant  settling  during  the  rise  and  help  simulate 
a  stabilized  cloud  stem. 

DELFIC  with  Vortex.  The  results  with  vortex  flow  modeling  show  a  drastically 
different  stabilized  cloud.  There  is  less  settling  for  large  particles  until  the  cloud  has  verti 
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cally  stabilized  and  the  toroidal  flow  is  turned  off.  The  majority  of  the  settling  for  the  vor¬ 
tex  modeling  occurs  after  the  cloud  has  vertically  stabilized.  This  does  allow  for  the  larger 
particles  to  maintain  a  much  higher  distribution,  with  virtually  nothing  having  settled  to 
the  surface  by  stabilization  time.  The  smaller  particles,  on  the  other  hand,  do  not  stabilize 
as  high  in  the  cloud  as  they  did  in  the  1-D  modeling.  This  too  is  due  to  the  vortex  flow 
moving  the  particles  about  in  a  toroidal  motion  instead  of  just  an  updraft  and  gravity  set¬ 
tling  as  in  the  current  model. 

Radial  Distribution  and  Fallout  Contours.  To  get  an  idea  of  the  effect  of  the  vortex 
flow  field  modeling  on  fallout  contours,  dose  rate  contours  were  generated  and  compared 
to  the  current  modeling.  The  current  DTM  and  OPM  employ  modeling  which  basically 
loses  the  radial  definition  of  the  2-D  stabilized  cloud.  An  empirical  radial  distribution  is 
superimposed  on  the  grounded  parcels.  The  difference  in  fallout  contours,  therefore,  only 
show  the  effect  of  the  change  in  the  vertical  distribution. 

The  results  of  this  type  of  comparison  for  the  1979  test  case  are  shown  in  Figure  42  on 
page  115.  The  one-hour  dose  rates  are  given  for  the  fallout  generated  from  both  the  1-D 
stabilized  cloud  and  the  vortex  flow  field  cloud.  The  main  difference  noted  is  that  the  100 
R/hr  contour  is  extended  over  10  kilometers  farther  downwind  with  the  vortex  stabilized 
cloud.  Also,  the  hot  spots  for  dose  rates  above  1000  R/hr  are  virtually  eliminated.  Both  of 
these  differences  are  due  to  the  fact  that  the  radioactivity  associated  with  the  larger  parti¬ 
cles  is  stabilizing  at  a  higher  altitude  and  therefore  is  spread  over  a  larger  area.  This  results 
in  increased  diffusion  from  the  wind  and  a  larger  radius  for  the  disks.  A  disk  which  starts 
transport  with  the  main  cloud’s  radius,  as  opposed  to  a  disk  which  starts  transport  below 
the  cloud  with  a  stem  radius,  has  its  radioactivity  spread  over  a  larger  surface  area. 
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Figure  42:  Fallout  Contour  Comparison  for  1979  Test  Case  (contours  represent  1, 10, 
100,  &  1000  R/hr  Dose  Rates  at  H+l  hr;  distances  are  in  meters) 

115 


lUlllllWiM 


This  chapter  presents  a  summary  of  the  findings  of  the  research  into  nuclear  cloud  rise 
and  growth.  The  chapter  contains  sections  for  bubble  rise  simulation,  particle  rise  simula¬ 
tion,  and  recommendations  for  future  research  into  these  areas. 

Bubble  Rise  Simulations 

The  main  conclusions  on  bubble  rise  concerns  the  performance  of  a  single  term 
entrainment  equation  and  the  use  of  constants  for  parameters  which  were  previously  yield 
dependent.  It  was  shown  that  the  single  term  entrainment  equation  (See  Eqn  (12))  could 
be  used  in  a  set  of  eight  nonlinear  ordinary  differential  equations  to  achieve  an  improved 
level  of  performance  in  predicting  the  heights  of  nuclear  cloud  tops.  The  values  for  the 
entrainment  parameter  and  the  eddy  viscous  drag  parameter  were  changed  to  constant  val¬ 
ues  in  the  prc acting  of  nuclear  cloud  rise  and  growth. 

The  original  CRM  of  DELFIC  contained  the  single  term  form  of  the  entrainment  equa¬ 
tion.  It  was  based  on  the  theory  presented  in  the  literature  then,  and  still  being  used  today. 
The  revision  in  1970  to  a  multiple  term  entrainment  equation  was  based  partly  on  model¬ 
ing  choices  of  the  cloud  shape.  This  change,  and  the  change  to  yield  dependent  parame¬ 
ters,  were  attempts  to  match  the  observation  of  nuclear  cloud  tops  by  abandoning  the 
accepted  form  of  the  equations.  When  the  modeling  choices  used  in  the  1970  derivation 
were  abandoned  in  1977,  the  entrainment  equation  should  have  been  returned  to  its  single 
term  form. 

The  entrainment  and  eddy  viscous  drag  parameters  are  the  two  most  important  model 
parameters  in  the  cloud  rise  equations.  Each  originated  as  constants  in  the  1967  version  of 


116 


the  CRM.  Both,  however,  have  been  changed  over  the  years  to  yield  dependent  parame¬ 
ters.  This  research  confirmed  that  the  increase  with  yield  in  the  entrainment  parameter  was 
being  offset  by  the  decrease  in  the  eddy  viscous  drag  parameter.  The  suggested  constant 
value  for  the  entrainmtu,  parameter  is  0. 12  and  the  value  for  the  eddy  viscous  drag  param¬ 
eter  should  be  0.10.  These  values  are  based  on  a  minimization  of  the  FRMS. 

Both  of  the  above  conclusions  were  drawn  after  a  careful  validation  using  actual  US 
atmospheric  nuclear  test  observations.  A  set  of  54  events  with  a  wide  range  of  yields  and 
heights  of  bursts  was  used.  These  computer  simulations,  using  recorded  atmospheric 
inputs,  showed  that  the  improvements  slightly  increase  the  ability  to  match  cloud  tops. 

As  a  preliminary  step  to  making  the  above  changes,  corrections  were  made  in  the 
CRM  subroutines.  Inconsistencies  with  the  stated  theory  in  the  1979  documentation  were 
removed.  These  corrections  also  made  a  slight  increase  in  the  performance  of  the  model. 

Eamclc-Rific  Simulations 

The  major  conclusion  drawn  in  the  particle  rise  studies  was  that  better  agreement  with 
hydrocodes  in  simulating  the  stabilized  radioactive  cloud  can  be  achieved  by  using  a  2-D 
flow  field.  The  particle  rise  modeling  of  the  CRM  in  1-D  was  replaced  with  numerical 
simulation  of  an  analytic  flow  field  based  on  Hill’s  spherical  vortex.  As  a  result,  the  larger 
particles  (>200  microns)  in  the  fallout  size  distribution  stabilize  higher  in  the  atmosphere. 
This  is  consistent  with  comparisons  to  full  hydrocode  simulations. 

1116  new  model  accounts  for  circulation  in  a  toroidal  field  within  the  rising  bubble. 
This  was  accomplished  by  applying  a  known  analytical  solution  of  a  rising  bubble  in  a  still 
atmosphere  to  nuclear  cloud  rise.  This  added  information  of  the  flow  field  was  achieved 
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without  resorting  to  a  full  hydrocode  solution.  The  characteristics  of  the  rising  bubble 
needed  for  the  solution  of  the  flow  field  were  provided  by  the  bubble  rise  simulations  dis¬ 
cussed  in  the  previous  section. 

By  applying  the  above  particle  rise  method  to  nuclear  clouds,  increased  fidelity  was 
achieved  in  determining  the  locations  of  the  radioactive  dust  in  the  stabilized  nuclear 
cloud.  Previous  studies  showed  that  there  was  an  inconsistency  in  the  stabilized  heights  of 
larger  fallout  particles  when  compared  to  the  dust  loading  produced  by  hydrocodes.  The 
reason  for  the  discrepancy  was  shown  to  be  in  the  way  the  CRM  simulates  the  panicle  rise 
in  one  dimension.  By  accounting  for  a  spatially  varying  velocity  field  which  includes  cir¬ 
culation,  these  larger  particles  do  achieve  a  higher  altitude.  The  circulation  in  the  cloud 
keeps  the  particles  within  the  cloud  longer,  thereby  allowing  them  to  achieve  a  higher  alti¬ 
tude. 

Comparisons  to  hydrocode  dust  loading  plots  indicate  that  the  radioactive  debris  in  the 
nuclear  cloud  concentrates  in  a  toroidal  core  within  the  rising  bubble.  This  spatial  parti¬ 
tioning  of  the  debris  which  starts  in  the  fireball  from  the  dust  which  enters  the  cloud  at 
later  times  is  seen.  When  comparing  dust  density  contours,  the  sweepup  dust  appears  to 
form  the  bulk  of  the  stem  and  fill  in  the  void  along  the  vertical  axis.  The  dust  that  is  swept 
up  into  the  cloud  at  later  times  may  not  be  as  radioactive  as  the  debris  in  the  torus.  The  tor¬ 
oidal  positioning  of  the  radioactive  dust  in  the  new  model  is  much  different  than  the 
bivariate  gaussian  distribution  assumed  by  the  DTM  and  OPM  modules  of  DELFIC. 

Changes  to  the  dose  rate  map  from  the  change  in  vertical  distribution  of  the  radioactiv¬ 
ity  were  noticed.  The  higher  altitudes  achieved  by  the  larger  particles  help  spread  the  asso¬ 
ciated  activity  over  a  larger  area.  This  reduces  the  number  of  hot  spots  from  the  early 
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downing  of  the  large  particles  but  increases  the  range  of  the  lower  intensity  contours.  The 

difference  in  radial  distribution  may  also  alter  the  dose  rate  contours. 

Recommendations 

1 .  As  hydrocodes  become  more  useful  and  their  output  more  trusted,  they  could  be  used 
to  initialize  the  soil  loading  in  the  particle  rise  of  DELFIC.  A  spatially  varying 
description  in  the  early  seconds  after  detonation  could  replace  the  assumed  uniform 
cylindrical  distribution  from  cloud  top  to  the  surface. 

2.  A  time  dependent  loading  of  the  cloud  should  be  investigated.  The  location  of  the  vor¬ 
tex  sphere  at  the  time  of  soil  loading  is  of  utmost  importance  in  determining  if  soil 
enters  the  toroidal  core.  Soil  entering  the  cloud  after  the  initial  time  may  not  be  highly 
radioactive.  It  should  be  flagged  and  only  passed  onto  the  transport  phase  for  dose  rate 
calculations  if  it  becomes  collocated  with  the  initial  radioactive  debris  and  thus 
acquire  radioactivity  by  condensation  or  agglomeration. 

3.  The  spheroidal  vortex  theory  should  be  investigated  to  perfectly  align  the  vortex 
region  with  the  visible  oblate  bubble.  This  may  raise  the  location  of  the  stabilized  dust 
even  further  by  raising  the  vortex  midplane  to  the  center  of  the  visible  bubble. 

4.  A  method  for  accounting  for  residual  circulation  in  the  vertically  stabilized  cloud 
should  be  added  to  better  match  reality.  This  may  be  possible  by  relating  the  circula¬ 
tion  to  the  turbulent  kinetic  energy  density,  or  a  characteristic  velocity,  which  are  cur¬ 
rently  calculated  in  the  bubble  rise  equations.  Also,  a  method  that  allows  for  diffusive 
spreading  of  the  particles  back  towards  the  void  near  the  axis  should  be  investigated. 

5.  A  review  of  the  cloud  shape  in  the  bubble  rise  equations  should  be  done.  A  physical 
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mechanism  for  including  the  independent  tracking  of  the  cloud  top  and  base  after  the 
cloud  top  has  reached  apogee  is  in  order.  The  inclusion  of  oscillation  for  large  yield 
clouds  and  a  late  time  collapse  would  allow  for  a  closer  modeling  of  real  nuclear 
clouds.  Allowance  for  a  changing  cloud  eccentricity  should  be  made. 

6.  A  more  comprehensive  validation  of  the  bubble  rise  model  should  be  done  once  a 
change  to  the  cloud  shape  modeling  has  occurred.  An  increase  in  the  number  observed 
events  can  be  accomplished  now  with  the  declassification  of  test  data  not  available  to 
researchers  earlier.  The  gap  in  the  yield  range  from  1 10  kt  to  3.5  MT  should  be  filled  in 
along  with  more  shots  from  the  Pacific  Proving  Grounds. 

7.  The  transport  and  output  processor  modules  of  DELFIC  should  be  reviewed  for  inclu¬ 
sion  of  a  2-D  description  of  the  stabilized  radioactive  dust  cloud.  Once  a  2-D  cloud  is 
accepted  by  a  new  DTM/OPM,  the  effects  of  both  the  revised  vertical  distribution  and 
possible  effects  from  the  radial  distribution  on  fallout  contours  may  be  found.  A  com¬ 
parison  to  actual  recorded  fallout  footprints  from  the  US  atmospheric  nuclear  tests 
would  be  in  order. 
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Appcadix  A;  Cloud  Togs 


The  cloud  tops  of  the  Single  Burst  Nuclear  Cloud  Database  (SBNCD)  91  Series  of 
simulations  are  discussed  in  this  appendix.  In  particular,  the  cloud  tops  computed  in  the 
MAZ/TASS  simulations  are  compared  to  the  cloud  tops  computed  with  the  Cloud  Rise 
Module  (CRM)  of  DELFIC.  The  values  used  for  the  MAZ/TASS  hydrocode  simulations 
arc  the  10*8  g/cc  (i.e.  approximate  visible  density)  DST+PBL  soil  mass  contours  (9).  Two 
different  values  are  used  for  th~  CRM  simulations:  bubble  rise  visible  cloud  tops,  and  par¬ 
ticle  rise  10'8  g/cc  soil  mass  contours. 

All  cases  were  run  with  the  same  atmospheric  sounding  used  in  the  MAZ/TASS  simu¬ 
lations.  These  atmospheric  properties  included  temperature,  pressure,  humidity,  and  winds 
for  various  altitudes.  The  comparisons  made  in  this  study  use  only  soil  mass  for  the  den¬ 
sity  contours.  It  must  be  remembered  that  the  addition  of  condensed  water  mass  to  the 

density  contours  would  only  increase  the  predicted  cloud  top  values  by  raising  the  lO"8  g / 
cc  contour.  The  bubble  rise  cloud  top,  on  the  other  hand,  would  not  change,  as  it  is  based 
on  the  rise  of  the  entire  visible  bubble. 

Shown  in  Table  12  on  page  128  are  the  computed  cloud  tops  for  the  9  near  surface 
bursts  compared  in  Chapter  IV.  The  comparison  shows  that  the  MAZ/TASS  cloud  tops 
overpredict  the  CRM  bubble  rise  values  for  7  of  the  9  cases.  The  difference  between 
MAZ/TASS  and  the  CRM  increases  with  SHOB  for  those  yields  with  multiple  SHOBs. 

The  CRM  particle  rise  results  show  that  the  visible  contour  is  always  below  the  bubble 
rise  cloud  top  value.  This  is  due  to  the  modeling  choice  of  matching  the  visible  bubble 
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cloud  top  to  the  vortex  sphere  top.  This  ensures  that  all  particles  in  the  vortex  region  will 
be  below  the  validated  bubble  rise  cloud  top  value.  This  also  implies  that  even  if  the  con¬ 
densed  water  contribution  to  density  contours  were  added  to  the  soil  density,  the  visible 
contour  would  not  exceed  the  bubble  rise  cloud  top. 

Since  the  bubble  rise  cloud  top  predictions  are  validated  against  US  atmospheric 
nuclear  test  observations,  the  conclusion  is  that  MAZyTASS  cloud  top  predictions  are  too 
high.  This  is  especially  true  for  larger  SHOBs. 


Table  12:  MAZ/TASS  vs.  DELFIC  Cloud  Tops 


Case 

Yield  (kt) 

SHOB 

(ft/ktl/3) 

TASS 
Cloud  Top 
(km) 

Bubble  Rise 
Cloud  Top 
(km) 

Particle 
Rise  Cloud 
Top  (km) 

91H01 

20 

0 

11.5 

8.0 

6.5 

91B01 

100 

0 

12.5 

13.5 

11.0 

91B02 

100 

50 

15.5 

13.7 

11.0 

91B03 

100 

120 

17.5 

13.8 

105 

91A01 

300 

0 

17.5 

15.6 

13.0 

91A02 

300 

50 

17.8 

15.7 

14.0 

91A03 

300 

120 

20.3 

15.9 

12.5 

91F01 

800 

0 

18.0 

18.0 

15.5 

91F02 

800 

50 

23.0 

18.3 

16.5 
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Appendix  B;  Latent  Heat 


The  energy  in  the  buoyantly  rising  bubble  is  affected  by  the  release  of  latent  heat  dur¬ 
ing  its  rise  to  lower  temperatures  and  pressures.  The  amount  of  latent  heat  released  in  the 
DELFIC  Cloud  Rise  Module  (CRM)  simulations  of  the  Single  Burst  Nuclear  Cloud  Data¬ 
base  91  Series  cases  is  discussed.  The  values  are  compared  to  the  much  larger  values  of 
the  MAZ/TASS  hydrocode  simulations. 

To  begin,  the  cloud  gets  its  initial  buoyant  energy  from  the  nuclear  detonation.  In  the 
CRM  of  DELFIC  this  is  modeled  as  45%  of  the  total  device  yield.  In  MAZ/TASS  32.5% 
of  the  yield  is  used  initially  (9).  In  addition,  as  the  cloud  cools  to  condensation  tempera¬ 
ture  for  the  vapor  in  the  cloud  an  additional  amount  of  energy  is  released.  This  latent  heat 
is  determined  in  the  CRM  simulations  by  the  ratio  of  condensed  water  in  the  cloud  to  that 
of  air.  By  determining  the  amount  of  air,  the  mass  of  condensed  water  in  the  cloud  is  deter¬ 
mined.  The  latent  heat  released  to  produce  this  amount  of  condensed  water  is  then  found. 

Table  13  on  page  130  shows  the  latent  heat  values  calculated  with  the  two  codes  on  the 
91  Series  cases.  The  values  listed  for  the  CRM  are  the  maximum  value  seen  during  the 
cloud  rise  and  growth  phase  of  development.  The  values  for  TASS  were  provided  by  Dunn 
(9).  All  simulations  used  the  same  atmospheric  sounding  as  input. 

The  CRM  values  indicate  a  trend  of  increasing  latent  heat  released  with  increasing 
yield  and  decreasing  SHOB.  The  TASS  values  do  not  show  such  a  trend.  The  TASS  simu¬ 
lations  also  show  a  much  larger  release  of  latent  heat  by  a  factor  from  almost  3  to  over  10. 
This  larger  latent  heat  release  is  somewhat  offset  in  the  total  effective  buoyant  energy  by 
the  CRM’s  use  of  a  larger  fraction  of  the  initial  device  yield.  However,  the  TASS  simula- 


129 


tions  still  calculate  a  much  larger  effective  buoyant  energy  for  all  cases  by  up  to  50%. 

The  reason  for  this  difference  may  lie  in  the  basic  assumption  of  the  single  cell 
approach  to  DELFIC  bubble  rise  modeling.  By  allowing  only  a  single  temperature  for  the 
whole  cloud,  the  entire  cloud  will  cool  and  condense  at  the  same  time.  In  a  hydrocode 
approach  to  the  solution,  the  cloud  is  allowed  to  keep  certain  cells  in  the  multicell  cloud 
hot  while  a  majority  of  the  cells  can  cool  to  much  lower  temperatures.  In  addition  to 
allowing  for  hot  spots  in  the  cloud,  the  bydrocode  solution  allows  for  local  concentrations 
of  water  vapor  and  therefore  condensation.  These  two  differences  may  allow  the  TASS 
simulations  to  condense  more  of  the  water  vapor  present  in  the  particular  atmosphere 
used. 


Table  13:  MAZ/TASS  vs.  CRM  Latent  Heat  and  Effective  Buoyant  Energy 


Case 

SHOB 

(ft./kt1/3) 

Yield  (kt) 

TASS 
Latent 
Heat  (kt) 

Tass 

Effective 

Buoyant 

Energy 

(kt) 

CRM 
Latent 
Heat  (kt) 

CRM 

Effective 

Buoyant 

Energy 

(kt) 

91H01 

0 

20 

18 

24 

7 

16 

91B01 

0 

100 

68 

100 

19 

64 

91B02 

50 

100 

50 

82 

18 

63 

91B03 

120 

100 

64 

97 

17 

62 

91B05 

400 

100 

42 

74 

13 

58 

91A01 

0 

300 

76 

174 

23 

158 

91A02 

50 

300 

120 

218 

21 

156 

91A03 

120 

300 

134 

231 

IS 

153 

91  AOS 

400 

300 

100 

197 

12 

147 

91F01 

0 

800 

175 

435 

25 

385 

91F02 

50 

800 

263 

523 

,  —  -  — — 

23 

. 

383 
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